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Abstract  A  continuum  mechanics  theory  of 
deformable  solids  is  formulated  to  account  for  large 
deformations,  nonlinear  elasticity,  inelastic  deforma¬ 
tion  mechanisms,  micro  structure  changes,  and  time 
dependent  fields,  i.e.,  dynamics.  The  theory  incorpo¬ 
rates  notions  from  Finsler  differential  geometry,  and 
it  provides  a  diffuse  interface  description  of  surfaces 
associated  with  microstructure.  Mechanisms  include 
phase  transitions  and  inelastic  shearing,  with  phase 
boundaries  and  shear  planes  the  associated  surfaces. 
A  director  or  internal  state  vector  of  pseudo-Finsler 
space  is  viewed  as  an  order  parameter.  Newly  derived 
in  the  present  work  are  the  governing  equations  for 
dynamics,  including  kinematic  relations,  balances  of 
momentum  and  energy,  and  evolution  law(s)  for  the 
internal  state.  Also  derived  are  jump  conditions  per¬ 
tinent  to  shock  loading.  Metric  tensors  and  volume 
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can  vary  isotropically  with  internal  state  via  a  confor¬ 
mal  transformation.  The  dynamic  theory  is  applied  to 
describe  shock  loading  of  ceramic  crystals  of  boron 
carbide,  accounting  for  inelastic  mechanisms  of  shear 
accommodation  and  densification  upon  amorphization 
under  high  pressure  loading.  Analytical  predictions 
incorporating  the  pseudo-Finsler  metric  demonstrate 
remarkable  agreement  with  experimental  data,  with¬ 
out  parameter  fitting.  Additional  solutions  suggest  that 
dynamic  shear  strength  could  be  improved  significantly 
in  boron-based  ceramics  by  increasing  surface  energy, 
decreasing  inelastic  shear  accommodation  in  softened 
amorphous  bands,  and  to  a  lesser  extent,  by  increasing 
the  energy  barrier  for  phase  transformation. 

Keywords  Continuum  mechanics  •  Differential 
geometry  •  Shock  physics  •  Dynamic  fracture  • 

Phase  field  •  Ceramics 


1  Introduction 

Crystalline  solids  may  demonstrate  a  number  of  inelas¬ 
tic  deformation  and  failure  mechanisms  under  intense 
dynamic  loading.  These  mechanisms  include  crystallo¬ 
graphic  slip  (i.e.,  dislocation  glide),  deformation  twin¬ 
ning,  adiabatic  shear  localization,  phase  transforma¬ 
tions,  dynamic  cleavage  and  intergranular  fracture,  and 
nucleation,  growth,  and  coalescence  of  voids.  Efforts 
towards  modeling  the  physics  of  such  phenomena  in  the 
context  of  continuum  constitutive  theory  and  numeri- 
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cal  simulations  have  witnessed  incremental  progress 
for  several  decades.  Models  for  continuum  crystal 
plasticity  at  high  strain  rates  include  (Clayton  2005a; 
Luscher  et  al.  2013;  Lloyd  et  al.  2014a;  Shahba  and 
Ghosh  2016).  Models  for  deformation  twinning  include 
pseudo- slip  approaches  (Barton  et  al.  2009;  Clayton 
2009)  and  phase  field  approaches  (Clayton  and  Knap 
2011a,  2013,  2015a;  Hildebrand  and  Miehe  2012). 
Theoretical  and  computational  continuum  mechanics 
models  have  been  developed  accounting  for  shear  local¬ 
ization  (Li  et  al.  2001;  Bronkhorst  et  al.  2006;  Sun  and 
Mota  2014)  and  ductile  damage,  i.e.,  void  mechanics 
(Bammann  and  Solanki  2010;  Bronkhorst  et  al.  2016). 
Phase  field  theory  has  also  been  applied  to  describe 
solid- solid  phase  transformations  in  the  context  of  finite 
deformations  (Levitas  et  al.  2009;  Clayton  2014a;  Lev- 
itas  2014).  Fractures  have  been  modeled  via  a  num¬ 
ber  of  computational  schemes,  including  cohesive  zone 
methods  (Xu  and  Needleman  1994;  Li  and  Wang  2004; 
Clayton  2005b;  Vogler  and  Clayton  2008;  Foulk  and 
Vogler  2010)  and  phase  field  representations  (Borden 
et  al.  2012;  Clayton  and  Knap  2014,  2015b). 

The  current  paper  invokes  a  Finsler-geometric 
theory  of  continuum  mechanics  (Clayton  2016a,  b, 
2017a,  b)  to  account  for  nonlinear  elasticity,  inelastic¬ 
ity,  and  microstructure  evolution  under  dynamic  load¬ 
ing  conditions.  In  the  context  of  geometrically  nonlin¬ 
ear  continuum  physics,  the  body  manifold  in  its  refer¬ 
ential  and  spatial  configurations  is  treated  as  a  Finsler 
space  rather  than  a  Euclidean  space  of  conventional 
continuum  mechanics.  The  metric  tensor  and  its  deriva¬ 
tive  quantities  such  as  lengths,  areas,  and  volumes  all 
functionally  depend  on  an  internal  state  vector  or  direc¬ 
tor  vector.  More  precisely,  in  the  present  class  of  theory 
(Clayton  201 6a,  b,  2017a,  b),  each  configuration  is  con¬ 
sidered  of  pseudo-Finsler  geometry  since  the  internal 
state  vector  need  not  be  a  unit  vector  and  the  metric 
need  not  be  homogeneous  of  degree  zero  with  respect 
to  internal  state  (Minguzzi  2014).  In  this  work,  as  first 
proposed  in  Clayton  (2016a,  2017a),  a  conformal  trans¬ 
formation,  i.e.,  Weyl  scaling  (Weyl  1952),  is  used  to 
account  for  dependence  of  the  metric  on  internal  state. 
Prior  work  in  the  context  of  quasi-static  loading  condi¬ 
tions  has  demonstrated  similarities  between  this  class 
of  Finsler  theory  (Clayton  2016a)  and  variational  phase 
field  models  (Clayton  and  Knap  2015a)  when  the  inter¬ 
nal  state  is  regarded  as  an  order  parameter.  Historically, 
work  in  Finsler  geometry  applied  to  solid  mechan¬ 
ics  is  limited  in  scope:  the  ferromagnetic  theory  of 


Amari  (1962),  kinematic  models  of  Bejancu  (1990), 
and  a  more  comprehensive  thermomechanical  theory 
with  a  numerical  example  in  Saczuk  (1996),  Stumpf 
and  Saczuk  (2000).  A  recent  review  of  the  subject  is 
Clayton  (2015a). 

As  explained  in  Clayton  (2017b),  generalized  Finsler 
geometric  continuum  mechanics  includes  some  dis¬ 
tinctive  features  absent  in  continuum  mechanics  mod¬ 
els  of  standard  Riemannian  geometry.  Pseudo-Finsler 
space  invokes  basis  vectors,  metric  tensors,  and  connec¬ 
tion  coefficients  that  may  explicitly  depend  on  the  inter¬ 
nal  state  or  director  vector.  Volume  and  area  elements 
inherit  an  explicit  state  dependence,  and  Stokes’  the¬ 
orem  includes  terms  accounting  for  state  dependence 
of  connection  coefficients  in  covariant  differentiation. 
Resulting  equilibrium  equations  thereby  contain  cor¬ 
respondingly  novel  terms  that  arise  from  application 
of  Stokes’  theorem.  This  enriching  of  the  governing 
equations  is  thus  motivated  by  mathematical  physics 
rather  than  ad-hoc  constitutive  assumptions.  Physi¬ 
cally,  dependence  of  the  metric  tensor  on  internal  state 
is  realistic  for  solid  bodies  when  evolving  microstruc¬ 
ture  alters  lengths  and/or  orientations  of  material 
line  elements.  Finsler  geometry  is  not  required  for 
mathematically  rigorous  modeling  of  micro  structure 
changes  or  inelastic  deformation,  but  it  does  provide  an 
enhanced  description  relative  to  standard  Riemannian- 
geometric  representations.  The  present  general  frame¬ 
work  also  has  been  shown  (Clayton  2017a)  to  encom¬ 
pass  and  reduce  to  (under  certain  simplifying  assump¬ 
tions)  other  classes  of  continuum  models  such  as  those 
of  micropolar  (Clayton  et  al.  2006),  strain  gradient 
(Clayton  et  al.  2004a),  and  phase  field  type.  Further 
advantages  may  become  evident  upon  examination  of 
solutions  to  specific  problems. 

The  present  paper  extends  previous  variational- 
based  (i.e.,  quasi-static  or  incremental)  theory  (Clay¬ 
ton  2016a,  b,  2017a,  b)  to  the  dynamic  regime.  Govern¬ 
ing  equations  for  dynamics  are  newly  derived:  balances 
of  mass,  linear  momentum,  and  energy  and  evolution 
equation(s)  for  the  internal  state.  The  latter  equation(s) 
are  similar  in  form  to  those  of  Ginzburg-Landau  or 
Allen-Cahn  phase  field  theory  (Allen  and  Cahn  1979; 
Levitas  2014),  albeit  with  additional  terms  manifest¬ 
ing  to  account  for  Finsler  space.  An  important  mathe¬ 
matical  device  used  in  the  current  derivations  centers 
on  the  divergence  theorem  of  Finsler  geometry  first 
presented  by  Rund  (1975).  Also  newly  addressed  is  a 
Finsler-geometric  description  of  shock  wave  propaga- 
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tion  in  solids.  Jump  conditions  are  derived  by  extend¬ 
ing  methods  in  Casey  (2011)  and  Clayton  (2014b,  c) 
for  Euclidean  space  to  pseudo-Finsler  space.  In  this 
context,  discontinuities  in  stress,  deformation  gradi¬ 
ent,  and  particle  velocity  are  assumed  to  exist  across  a 
shock  front  moving  at  steady  velocity.  A  linear  differ¬ 
ence  approximation  is  used  to  represent  state  variable 
gradients  across  the  front,  and  regions  far  behind  the 
front  are  assumed  to  be  in  equilibrium  with  respect  to 
internal  state. 

The  Finsler  description  of  shock  physics  is  used  to 
model  planar  impact  in  boron  carbide  crystals.  Boron 
carbide  (B4C)  is  a  ceramic  material  of  keen  current 
interest.  Beneficial  physical  properties  for  industrial 
and  defense  applications  include  high  hardness,  high 
elastic  stiffness,  and  low  mass  density.  Detrimental 
properties  include  low  ductility  and,  under  impact  load¬ 
ing  conditions,  a  tendency  for  failure  by  shear  localiza¬ 
tion.  Stress  induced  amorphization  is  thought  to  pre¬ 
cede  and/or  accompany  this  failure  mechanism.  Shear 
banding  is  a  prominent  inelastic  deformation  mecha¬ 
nism  in  metallic  glasses  as  well  (Cheng  and  Ghosh 
2013).  The  mechanism  of  amorphization  in  B4C  is 
thought  related  to  collapse  of  the  rhombohedral  unit  cell 
under  compression  along  the  c-axis  ([0001]  in  hexag¬ 
onal  Miller  indices)  and  is  further  enabled  by  shearing 
modes  (Grady  201 1 ;  Taylor  et  al.  2012;  Clayton  2012a, 
2013a,  2014a;  Taylor  2015;  An  and  Goddard  2015a). 

Prior  work  applied  (Clayton  2016b,  2017b)  a  static 
version  of  Finsler  theory  to  describe  tension,  compres¬ 
sion,  and  shear  of  boron  carbide  single  crystals.  In  this 
paper,  the  dynamic  version  of  Finsler  theory  is  newly 
applied  to  shock  compression  of  this  ceramic.  An  order 
parameter  is  linked  simultaneously  to  densification  and 
inelastic  shear  accommodation  as  the  crystal  transitions 
to  a  glassy  phase  under  intense  dynamic  loading.  Ana¬ 
lytical  solutions  are  derived  in  terms  of  a  set  of  algebraic 
equations  that  may  be  solved  simultaneously  (albeit, 
not  in  closed  form)  for  the  shock  stress,  order  param¬ 
eter,  entropy,  and  shock  characteristics  such  as  shock 
velocity  and  particle  velocity.  In  addition  to  model¬ 
ing  the  response  of  the  nominal  material,  the  current 
solutions  are  used  to  probe  effects  of  changes  in  fun¬ 
damental  material  properties  on  dynamic  shock  stress 
and  dynamic  shear  strength,  the  latter  of  high  impor¬ 
tance  in  the  context  of  resistance  of  ceramic  materials  to 
ballistic  penetration  (Curran  et  al.  1993;  Bourne  2008; 
Clayton  and  Tonge  2015;  Clayton  2015b,  2016c).  Such 
studies  are  motivated  by  concurrent  efforts  at  computa¬ 


tional  materials  design  (e.g.,  via  atomic  simulations)  of 
boron-based  ceramics  for  improved  ductility  and  fail¬ 
ure  resistance,  for  example  alloying  with  boron  subox¬ 
ide  (An  and  Goddard  2015b;  Tang  et  al.  2015). 

The  remainder  of  this  paper  is  organized  as  fol¬ 
lows.  Section  2  presents  the  general  theory  for  contin¬ 
uum  dynamics  in  the  absence  of  discontinuities,  includ¬ 
ing  differential  geometry,  kinematics,  thermodynam¬ 
ics,  and  balance  laws.  Section  3  develops  the  conser¬ 
vation  or  jump  conditions  for  a  solid  domain  through 
which  a  shock  wave  propagates  at  a  steady  material 
velocity.  Section  4  reports  application  of  the  dynamic 
theory  to  planar  shock  loading  of  B4C.  Section  5  con¬ 
tains  concluding  remarks. 

Notation  of  continuum  physics  and  Finsler  geome¬ 
try  is  used.  For  coordinate-free  descriptions,  bold  type 
is  used  for  vectors  and  tensors,  and  italic  type  is  used 
for  scalars  and  scalar  components.  Index  notation  is  fre¬ 
quently  invoked  with  the  usual  Einstein  convention  for 
summation  over  repeated  contravariant  and  covariant 
pairs. 

2  Theory 

Finsler-geometric  descriptions  of  reference  and  spatial 
configurations  are  given,  followed  by  a  description  of 
motion  of  the  body.  Mathematical  aspects  discussed 
already  in  prior  work  (Clayton  2016b,  2017  a,  b)  are 
presented  only  in  brief.  Kinematics  are  further  devel¬ 
oped  to  account  for  time  dependency,  e.g.,  velocities 
and  deformation  rates.  Thermodynamic  principles  and 
corresponding  balance  laws  are  then  newly  defined  or 
derived  for  the  dynamic  case. 

2.1  Pseudo-Finsler  geometry 

Finsler  geometric  representations  of  reference  and  spa¬ 
tial  configurations  are  given  in  Sects.  2.1.1  and  2.1.2. 
Mappings  between  such  configurations  are  then  dis¬ 
cussed  in  Sect.  2.1.3.  Rate  kinematics  are  defined  in 
Sect.  2. 1 .4.  Multiplicative  decompositions  of  the  defor¬ 
mation  gradient  and  metric  tensors  are  introduced  in 
Sect.  2.1.5. 

2.1.1  Reference  configuration 

Fet  201  be  a  differential  manifold  of  spatial  dimen¬ 
sion  3.  Fet  X  e  201  denote  a  material  point,  and  let 


<£)  Springer 


56 


J.  D.  Clayton 


{Xa}(A  =  1,  2,  3)  denote  a  coordinate  chart  covering 
part  or  all  of  9JT.  To  each  point  is  assigned  a  vector  D 
serving  as  a  descriptor  of  possibly  evolving  microstruc- 
ture.  Coordinates  {Da}(A  =  1,  2,  3)  are  the  entries  of 
D.  The  vector  Z>,  called  a  state  vector  or  internal  state 
vector,  need  not  be  of  unit  length,  thought  it  may  be 
referred  to  as  a  director  vector  regardless.  In  notational 
arguments,  dependence  of  a  function  on  (X ,  D)  implies 
dependence  on  reference  charts  ({XA},{DA}). 

Holonomic  basis  vectors  are  the  fields  {^J-4 ,  }. 

Let  D)  denote  nonlinear  connection  coeffi¬ 

cients.  In  Finsler  geometry,  the  non-holonomic  basis 
whose  entries  transform  between  coordinate  systems 
as  typical  vectors  is 


8 

8X A 


— J-na— R’  8Da  =dDA  +  N$dXB. 
dXA  AdDB  B 

(2.1) 


The  set  {  ,  ^4 }  is  logically  used  for  a  local  basis 

on  the  referential  tangent  bundle,  and  the  reciprocal 
set  {dXA,  8Da}  for  the  cotangent  bundle.  The  Sasaki 
metric  is 


G(X,  D)  =  Gab(X ,  D)dXA  0  dXB 

+  GAb(X ,  D)8Da  0  8Db.  (2.2) 


Components  Gab  (Gab)  are  used  to  lower  (raise) 
indices,  and  the  determinant  is  G(X,  D)  = 
det [Gab  (X,  D)] .  Partial  coordinate  differentiation  and 
delta-differentiation  are  written  as 


9a(  •) 


«a(0 


9(0  -  9(0 

'dXA'  9a</  ~  '()DA  ’ 

=  dA(.)  -  NBdB(-). 


(2.3) 


The  Christoffel  symbols  of  the  second  kind  for  the  Levi- 
Civita  connection  on  are 


Ybc  =  \GAD (8cGbd  +  $bGcd  ~  9 dGbc ) 

=  Gadybcd •  (2.4) 


Cartan’s  tensor  referred  to  reference  coordinates  is 


The  spray  and  its  canonical  nonlinear  connection  coef¬ 
ficients,  with  the  latter  an  example  of  those  in  (2.1) 
when  =  G^,  are,  respectively, 

Ga  =  \ybCDbDc  ,  =  dBGA.  (2.7) 


Denote  by  V  the  co variant  derivative.  Horizontal  gra¬ 
dients  of  basis  vectors  are 


V S/SXB  -  H 

9 

vs/sxB  JBC  -  K 


BC  8XA’ 
A  J_ 
BC  dDA’ 


(2.8) 


with  coefficients  H  AjC  and  KgC.  Vertical  gradients  of 
basis  vectors  are  similarly  obtained  from  coefficients 


VBC  an<^  ^BC: 

d 

T  T  A 

d 

V3/30B  qdc  ~ 

^BC 

dDA ’ 

8 

,  ,  A 

8 

va/3D^zc  - 

ybc 

8XA' 

(2.9) 


The  above  descriptions  pertain  to  both  pseudo-Finsler 
space  and  Finsler  space.  The  latter  type  of  space  is  a 
subset  of  the  former,  and  conditions  for  which  a  pseudo- 
Finsler  space  is  a  strict  Finsler  space  are  discussed  else¬ 
where  (Minguzzi  2014;  Clayton  2016a).  Two  connec¬ 
tions  often  encountered  in  the  (pseudo)-Finsler  litera¬ 
ture  are  (Minguzzi  2014) 


Chern-Rund  connection: 


K$c 


Na 


=  GA,HA 


yA  —  yA  _  n. 
V  r  —  1  n r  —  U, 


BC 


~  ^ BC  ’  y  BC  ~  ±BC 
Cartan  connection:  = 


A  tj  A  _  jr  A  _ 

n  •  1 1 B C  —  nr  — 


r \ 


yA  _  yA  _  yA 
v  nr  —  1  BC  ~  ^ 


BC 


BC ’  v  BC 


BC • 


Let  (*)|c  denote  horizontal  covariant  differentiation  in 
a  coordinate  chart  {Xc }.  Then  when  either  of  these  two 
connections  is  used,  the  horizontal  covariant  derivative 
Gab\c  =  0. 

Denote  by  dX  and  and  d D  differential  line  or  vector 
elements.  Squared  differential  line  lengths  with  respect 
to  (2.2)  are 


CAc  =  \GAD(dcGBD  +  ~dBGCD  ~  dDGBC ) 

=  GadCBcd •  (2.5) 


|dZ|2  =  (dZ,  GdZ>  =  Gab dXAdXB , 
\dD\2  =  (d D,  GdD)  =  GABdDAdDB. 


(2.10) 


The  horizontal  coefficients  of  the  Chern-Rund  and  Car- 
tan  connections  are  both  equal  to 

r 'bc  =  \GAD (8cGbd  +  8bGcd  -  &dGbc) 

=  GADrBCD.  (2.6) 


Scalar  volume  elements  and  volume  forms  of  9JI  are 
Rund  (1975) 

dv  =  v/GdZ1dZ2dZ3, 

d£2  =  VGdZ1  A  dZ2  A  dZ3.  (2.11) 
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The  differential  area  form  corresponding  to  a  compact 
region  of  3DT  is 

S2  =  /fidU1  Ad[/2;  [xA  =  XA(Ua),  (a  =  1,2); 


o/ 

Pa 


dxA 

3E/“’ 


P  =  det  (pAGABpJj) 


(2.12) 


The  following  identities  apply  for  gradients  of  G: 

SA(lnVG)  =  rbB, 

(VG)|A  =  9a(VG)  -  N*~dB(VG)  -  Vgh£b. 

(2.13) 

Let  VA(X,  D)Q(X,  D)  be  a  2-form,  and  let  VA  be 
contravariant  components  of  generic  differentiable  vec¬ 
tor  field  V  =  VA  .  Prescribe  the  horizontal  connec¬ 
tion  as  one  for  which  H%B  =  r^B  =>►  (VG)|a  =  0. 
Then  in  a  coordinate  chart  {XA},  Stokes’  theorem  is 
Rund  (1975) 


/ 

Jm 


[V,A  +  (V  Cbc  +  dBVA)D*A]dn 


vanaq. 


(2.14) 


am 


Here,  NA  is  the  unit  outward  normal  to  3  9JI,  VAA  = 
SaVa  +  VaH%a,  and  D°A  =  dADB  +  N%.  Motivated 
by  this  theorem,  a  covariant  derivative  operation  (*)||A 
is  defined  in  a  reference  coordinate  chart  as 

( ■ )  1 1  a  =  0)|  a  +  [(-)C£c  +  3b(-)]£>.Ba 


[  (Oil Ad£2  =  (f 
Jon  Jdo 


(0  NAC2.  (2.15) 

Jon  Jdon 

In  other  words,  (2.15)  enables  presentation  of  Stokes’ 
theorem  (including  the  divergence  theorem  of  Rund 
Rund  (1975))  in  Finsler  geometry  in  a  compact  form 
comparable  to  that  in  Riemannian  geometry  (Clayton 
2014d). 


2.1.2  Spatial  configuration 

The  spatial  configuration  corresponds  to  a  time  t  at 
which  the  solid  has  undergone  motion.  A  pseudo- 
Finsler  geometric  framework  exists  that  exactly  par¬ 
allels  Sect.  2.1.1.  Notation  provides  the  distinction: 
deformed  coordinates  and  their  indices  are  written  in 
lower-case  rather  than  capitals. 

Differential  manifold  m  of  spatial  dimension  3 
describes  the  body  embedded  in  ambient  Euclidean 
3 -space.  Let  v  e  m  denote  a  spatial  point,  and  let 
{xa}(a  =  1,  2,  3)  denote  a  chart  on  m.  Attached  to 


each  point  is  the  internal  state  vector  d ,  with  secondary 
coordinates  { da  }(a  =  1 ,  2,  3) ;  d  need  not  be  a  unit  vec¬ 
tor.  The  natural  or  holonomic  basis  is  {-^a,  jja  }•  With 
n%( x,  d)  spatial  nonlinear  connection  coefficients,  the 
usual  non-holonomic  basis  vectors  are 


8 

8xa 


8 

d¥’ 


8da  =  dda  +  nabdxb. 

(2.16) 


The  set  { ^ ^ }  will  be  used  as  a  local  basis  over  the 
spatial  tangent  bundle,  and  {dxa ,  8 da }  for  the  cotangent 
bundle.  The  spatial  Sasaki  metric  tensor  is 


g(x,d)  =  gab(x,d)dxa®dxb  +  gab(x,d)Sda®8db, 

(2.17) 


with  determinant  g(x,d)  =  det [gab(x,d)].  Spatial 
differentiation  follows  the  notation 


9a  (0 


9(0 

dxa  ’ 


da  (0 


9(0. 
3  da' 


8a(-)  =  j2  =  M-)~nbM-)-  (2-18) 

Levi-Civita  connection  coefficients  on  m  are 
Ybc  =  \gad(dcgbd  +  dbgcd  ~  ddgbc)  =  g^Ybcd- 


(2.19) 


Cartan’s  tensor  is 


Cbc  =  \gad(dcgbd  +  hgcd  ~  ddgbc)  =  gadCbcd. 

(2.20) 


Horizontal  coefficients  of  Chern-Rund  and  Cartan  con¬ 
nections  are 

rgc  =  \gad(tcgbd  +  hgcd  -  sdgbc )  =  gadrbcd. 

(2.21) 


The  spray  and  canonical  nonlinear  connection  coeffi¬ 
cients  (when  n%  =  gb)  are 

ga  =  \Ybcdbdc,  gab  =  3  bga.  (2.22) 


Horizontal  gradients  of  basis  vectors  are  found  from 
connection  coefficients  Hgc  and  K%c: 


'S/Sx1 


^-  =  K  — 


8xc 


8xa 


U/8x‘ 


—  Va 

3  d‘~Kbc 


3  da 
(2.23) 
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Vertical  connection  coefficients  Vfc  and  Ybc  yield  ver¬ 
tical  gradients: 


yd/ddb 


77  =  Vt 


ddc 


be 


dda 


/d/ddb 


8xc 


lbc 


8xa 

(2.24) 


Line  and  volume  elements  analogous  to  those  in 
Sect.  2.1.1  exist,  as  does  an  area  form  co  like  £2  in  (2.12). 
Spatial  versions  of  Stokes’  theorem  (2.14)  hold,  where 
in  particular,  the  analog  of  (2.15)  is 

(')||a  =  (0|a  +  [(•)  Ccbc  +  hi')  ]< 

=>  f  (Oil adco  =  ®  (0 na(0.  (2.25) 

2m  J  3m 


f(x,d,t )  =  fA(x,d,t) 


SXA 


(g)  dxa 


d<PA(x,d,t )  8 


dxa  8XA 
d  X(x,d,t)  A 


(g)  dx° 


dx 


fa  =  W*  =  daXA- 

(2.29) 


Functions  F  and  /  are  invertible  with  positive  deter¬ 
minants,  and  are  inverses  of  one  another  at  coincident 
points  on  VJl  or  m  at  the  same  time  instant  t: 


F%[X,  D(X,  0,  t]f£[ x(X,  D,  0,  d(x(X,  D,  t),  t)]  =  8ab , 
fA[x,  d( x,  0,  t]F%[X(x,  d,  0,  D(X(x ,  d,  t),  t)]  =  8$. 

(2.30) 


2.1.3  Deformation  kinematics 

The  time-dependent  motions  from  9Ji  to  m  and  vice 
versa  are  the  C2  functions 

xa(X,  D,t)  =  q>a[X,  D(X,t),t], 

XA(x ,  d,  t)  =  0A[x ,  d(x,  0,  f].  (2.26) 

Unlike  variational  theories  in  prior  work  (Clayton 
2016b,  2017a,  b),  here  time  t  is  an  explicit  independent 
variable.  Finsler  kinematics  (Bejancu  1990;  Saczuk 
1996;  Stumpf  and  Saczuk  2000;  Clayton  2017a)  may 
differ  from  classical  finite  kinematics  in  Riemannian 
geometry  via  incorporation  of  internal  state  in  these 
motion  functions.  State  vector  mappings  are 

da(X,D,t)  =  6a[X,D(X,t),t], 

DA(x,d,t)  =  0A[x,d(x,t),t].  (2.27) 

The  deformation  gradient  of  Finsler-geometric  con¬ 
tinuum  mechanics  (Clayton  2016b,  2017a)  is  defined 
as  the  partial  derivative  of  motion  referred  to  the  non- 
holonomic  basis: 


Director  deformation  gradients  introduced  elsewhere 
(Clayton  2016b)  are  admissible  but  inessential  to  forth¬ 
coming  derivations. 

Similar  to  those  of  Riemannian  geometry,  transfor¬ 
mation  equations  for  line  elements  and  volume  ele¬ 
ments/forms  are  (Clayton  2016a,  b,  2017a,  b) 

3x 

dx  =  — dX  dxa  =  Fa>  dXA, 
dX 

dX  =  —  dx  &  dXA  =  fAdxa\  (2.31) 

dx  a 

dv  =  JdV  =  [det(F^)7g7G]dV, 

dV  =  jdv  =  m(fA)^G/^]dv, 

dco  =  JdQ ,  d£2  =  jdco.  (2.32) 

Lengths  of  deformed  to  initial  line  elements  are 
described  by  the  symmetric  deformation  tensor  C  : 

Id*  I2  =  FaAFbBgabdXAdXB  =  CAB  dXAdXB 
=  (dZ,  CdX), 

C  =  CABdXA  ®dXB  =  FaAgabFbBdXA  ®  dXB . 

(2.33) 


F(X,  D,  t)  =  FUX,  D,  t) - <g>  dX 


Sxa 

d(pa(X,  D,  t)  8 

dX A  8x“ 
3  x(X,D,t) 

dX  ’ 


)dXA 

Fa  =  d  A<pa  =  dAXa. 

(2.28) 


The  inverse  tangent  mapping  from  spatial  to  referential 
coordinates  similarly  is 


It  follows  that  det(C^)  =  det (Cab)/G  =  J2. 
Mixed  horizontal  gradients  of  non-holonomic  bases  are 
obtained  using  (2.3),  (2.23),  and  (2.28): 

^sl*xAJxc  -  Jx^  Vs/Sx°  8^ 

£ 

=  (F“-NBdBxa)yS/Sxa  — 

=  {FaA-NBdBxa)Hbac-X.  (2.34) 
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2.1.4  Material  time  derivatives 


The  notion  of  a  material  time  derivative  in  Finsler  space 
is  needed  for  subsequent  derivations  in  the  dynamic 
Finsler-geometric  continuum  theory.  Let  a  superposed 
dot  and  the  notation  D(  )/Dt  denote  the  material  time 
derivative,  which  is  defined  here  as  the  partial  time 
derivative  of  a  quantity  at  a  fixed  material  point  X  and 
at  fixed  internal  state  D.  For  example,  let  A(X,  D,  t) 
denote  a  generic  differentiable  field  variable;  its  mate¬ 
rial  time  derivative  is 

DA(X,  D,  t)  dA(X ,  D,  t) 


A(X,  D,  t) 


Dt 


dt 


(2.35) 


The  material  velocity  (vector)  is  defined  as  the  material 
time  derivative  of  position,  i.e., 

dx(X,D,t)  d(pa(X,  D,  t)  8 


v(X,  D,t) 


dt 


dt 


Dxa 

~Dt 


8xa  ’ 
(2.36) 


The  material  acceleration  is  the  material  time  derivative 
of  velocity: 

dv(X,D,t)  n  8 

a(X,  D,  t)  =  ’  =  aa(X ,  D ,  t)  —  , 

dt  8xa 


Dva 

~Dt 


(2.37) 


A  fundamental  assumption  is  that  material  time  deriva¬ 
tives  of  Sasaki  metric  tensors  vanish  identically,  anal¬ 
ogously  to  the  same  identities  that  can  be  derived  in 
Riemannian  geometry  (Clayton  2014d)  where  metrics 
may  depend  only  on  coordinates  but  not  internal  state: 
dG(X,  D) 

G(X,  D)  =  - v  =  0, 

dt 

g[x(X,D,t),d(x(X,D,t),t)] 

_  dg[x(X,  D ,  t),  d(x(X,  D ,  t),  t)] 

~  dt 


0.  (2.38) 


This  description  naturally  excludes  convected  coordi¬ 
nate  representations.  The  first  of  (2.38)  is  derived  triv¬ 
ially  by  inspection;  the  second  assumes  that  a  metric- 
compatible  connection  describes  spatial  gradients  of 
basis  vectors  and  that  d  is  held  fixed  during  time  dif¬ 
ferentiation. 

The  definition  of  the  material  time  derivative  of 
the  deformation  gradient,  like  that  of  F  itself 
(Bejancu  1990;  Stumpf  and  Saczuk  2000;  Clayton 
2016a,  2017a),  is  non-unique  in  Finsler-geometric  con¬ 
tinuum  mechanics.  In  the  present  theory,  the  rate  of 
deformation  gradient  is  defined  as  follows: 


F(X,  D,  t)  = 


D 

~Dt 


F(X,  D ,  t )  =  L(X ,  D ,  t)F(X ,  D ,  t), 


(2.39) 


where  L  is  the  velocity  gradient  tensor: 


8  5 

L  =  Lci  - — j  (g)  dxb  =  Uml - 0  dxb . 


}8xa 


(2.40) 


The  covariant  derivative  operation  entering  definition 
(2.40)  is  defined  in  (2.25).  From  (2.38)  to  (2.40),  the 
time  derivative  of  the  Jacobian  determinant  of  (2.32)  is 


dJ 


J  =  ~aFl  =  J(F~YaFaA  =  JL“  =  Jv\ \a.  (2.41) 


Let  0  denote  the  volume  integral  of  a  generic  scalar 
quantity  0  over  spatial  domain  co ,  with  material  time 
derivative  4> : 


t)dco , 


<P(t)  =  [  (/>(x,D, 
J  O) 

d  r 

0(0  =  7)7  / 

F>t  J0J 


t)dco. 


(2.42) 


The  second  of  (2.42)  can  be  manipulated  using  (2.25) 
and  (2.41): 


D 

r 

0  =  - 

/  < 

Dt 

1  ^ 

JX2 

-i\ 

"30 

~L\ 

37 

bJdQ 


=  [(</> 

J  CD 


(0  +  0U|J  a)d(0 


\  x,D 


+  (0^fl)||a 


dco 


dco 


((j)va)naco. 


(2.43) 


In  this  version  of  Reynolds’  transport  theorem  extended 
to  Finsler  space,  the  partial  derivative  of  0(jc,  D,  t)  at 
fixed  x  and  D  is  defined  as  the  difference  30/3 1  =  0  — 
0||aiA  Equations  (2.35)-(2.43)  reduce  to  their  usual 
counterparts  of  classical  continuum  physics  (Clayton 
2011,  2014d)  when  D  and  d  are  omitted  entirely,  in 
which  case  the  respective  reference  and  current  config¬ 
uration  manifolds  become  those  of  usual  Riemannian 
geometry  and  Euclidean  space  in  particular. 


2.1.5  Multiplicative  kinematics 

A  multiplicative  decomposition  of  the 
Finsler-geometric  deformation  gradient  of  (2.28)  is 
invoked.  Let  F  be  decomposed  into  a  product  of  two 
non-singular  two-point  tensors: 
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F  =  FeFd,  FaA  =  ( FE)aa(FD)aA .  (2.44) 

Both  the  thermoelastic  or  thermomechanically  recov¬ 
erable  deformation,  FE,  and  the  residual  or  inelastic 
deformation  associated  with  changes  of  internal  state, 
Fd ,  have  positive  determinants.  Functional  forms  are 
(Clayton  2016b,  2017b) 

(. FEra  =  ( FEra(x ,  o,  ( FDrA  =  ( fd)ua[D(x ,  ol 

(2.45) 

where  the  inelastic  two-point  tensor  FD  may  have  fur- 
ther  dependence  on  X  only  via  possible  dependence  of 
its  basis  vectors  on  X.  Greek  indices  correspond  to  a 
generally  anholonomic  space  (Clayton  2012b,  2014d) 
(i.e.,  incompatible  intermediate  configuration).  Local 
integrability  conditions,  or  lack  thereof,  for  the  total 
deformation  gradient  and  its  constituents  are  discussed 
in  Clayton  (2016b).  A  multiplicative  split  of  the  director 
gradient  function  has  also  been  introduced  in  prior  work 
(Clayton  2016b).  As  demonstrated  in  detail  in  Clayton 
(2017b),  the  multiplicative  decomposition  of  deforma¬ 
tion  is  useful  for  partitioning  elastic  and  inelastic  contri¬ 
butions,  whereby  only  the  former  directly  influence  the 
mechanical  stress.  The  inelastic  deformation  may  gen¬ 
erally  consist  of  contributions  from  various  structural 
transformations  including  plasticity,  twinning,  phase 
changes,  and  damage  development.  Setting  FD  as  an 
explicit  function  of  D  enables  a  kinetic  relation  or  equi¬ 
librium  equation  for  the  internal  state  to  dictate  the  evo¬ 
lution  of  inelastic  deformation,  eliminating  the  need  to 
introduce  a  separate  governing  equation  for  the  latter. 

Multiplicative  decompositions  of  the  referential 
Sasaki  metric  tensor  and  its  inverse  are  used  later: 

Gab(X,D)  =  Gac(X)Gcb(D), 

Gab(X,  D)  =  Gac (X)(G_1)£(Z)).  (2.46) 

Analogous  decompositions  into  position-  and 
microstructure-dependent  parts  apply  for  the  spatial 
metric: 

8ab(x,d )  =  gac(x)gcb(d), 

gab(x,d)  =  gac(x)(g-l)c(d).  (2.47) 

Another  metric  tensor  is  used  on  the  intermediate  con¬ 
figuration  (Clayton  et  al.  2004b),  split  multiplicatively 
(Clayton  2016b): 


gap(X,D)  =  gay(X)gy(D), 

gaf>{X,D)  =  gay(X)Cg-X)PY(D)- 

g  =  det(ga^)  =  1/ det(g“^).  (2.48) 

For  the  total  intermediate  metric  and  intermediate 
structure-independent  metric,  the  following  forms 

(Clayton  2016b)  prove  most  convenient: 

gap(X,  D )  =  SaGab(X,  D)SB  =  SAGAB(X)gy(D)SB, 
gaft(X)  =  8AGAB(X)SB;  (2.49) 

A  local  volume  element  dv  and  volume  form  dco  on  the 
intermediate  space  are  obtained,  similarly  to  (2.32),  as 

dv  =  {det[(FD)aA]^/G}dV  =  JD  dV, 
d5  =  {det[(F£-1)“]y!7?}dv  =  jEdv; 
dd)  =  JD d£2  =  jEda).  (2.50) 

Jacobian  determinants  are  defined  for  inelastic  and  ther¬ 
moelastic  mappings: 

jD  =  Jd=  {det[(FD)«  ]7|7G)  = 

jE  =  Je=  {det[(JF£-1)“]7|7^}  =  (2-51) 

For  the  convenient  selection  (2.49),  g  =  G.  The  inelas¬ 
tic  deformation  FD  of  (2.44)  and  (2.45),  in  indicial 
notation  with  basis  vectors,  is 

Fd[D(X,  t)]  =  ( FDrB[D(X ,  t )]  ga  <g>  dXB 

=  ( FD)A[D(X,t)]gaAga®dXB , 

(2.52) 

with  ga  generic  basis  vectors  for  the  intermediate  state 
and  g shifter  components  from  reference  to  interme¬ 
diate  configurations  (Clayton  2014d).  For  consistency 
with  (2.49), 

8a=S«S^;  8a=Sa’  det(^)  =  1-  (2-53> 

2.2  Balance  equations  and  thermodynamics 

General  conservation  laws  for  Finsler-continuum 
dynamics  are  developed  in  Sect.  2.2.1.  Constitutive 
assumptions  leading  to  thermodynamic  identities  fol¬ 
low  in  Sect.  2.2.2.  Further  derivations  in  the  context  of 
the  multiplicative  kinematic  framework  of  Sect.  2.1.5 
conclude  the  section  in  Sect.  2.2.3. 
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2.2.1  General  conservation  principles 

Conservation  laws  for  mass,  momentum,  and  energy 
are  derived  in  what  follows  next.  This  treatment  extends 
the  quasi-static  variational  approach  of  prior  work 
(Clayton  2016b,  2017a,  b)  to  the  dynamic  regime  and 
reduces  to  the  latter  when  the  body  is  in  static  equilib¬ 
rium. 

Let  m  denote  the  mass  of  a  material  body,  and  dm  a 
differential  element  of  mass,  where 

d m(X)  =  po(X,  D)dV(X,  D) 

=  p[x(X ,  D,  t),  D ,  t] dv[x(X,  D ,  t),  D ,  t]. 

(2.54) 

Referential  mass  density  po  is  related  to  spatial  mass 
density  via  application  of  (2.32): 

p0  =  pj.  (2.55) 


a  local  angular  momentum  balance  (Clayton  2016b, 
2017a, b): 

aab  =  jgacPAFbA=oba.  (2.59) 

In  other  words,  Cauchy  stress  a  is  symmetric.  Let  B 
denote  a  local  body  force  vector  measured  per  unit  ref¬ 
erence  volume.  The  global  balance  of  linear  momentum 
is  posited  for  material  domain  3DT  as 

—  f  p0vdV  =  f  BdV  +  (f  t0dA.  (2.60) 
Dt  Jmt  Jmt  Jdvn 

Expressing  the  reference  traction  in  terms  of  Piola- 
Kirchhoff  stress  via  (2.58),  using  Stokes’  theorem 
(2.14)  for  generalized  pseudo-Finsler  space,  and  local¬ 
izing  the  result  to  a  differential  volume  element  at  point 
X  gives  the  following  local  balance  of  linear  momen¬ 
tum: 


As  in  classical  continuum  mechanics  in  the  absence 
of  mass  transport,  po  is  presumed  constant  in  time  at 
fixed  point  X ,  though  it  may  generally  vary  with  D. 
Therefore,  its  material  time  derivative  vanishes,  leading 
to  a  local  mass  balance  upon  use  of  (2.41): 

2-po  =  0  =►  P  +  pv“v,  =  yt+  \a  =  0.  (2.56) 

The  rightmost  equality  can  alternatively  be  derived  by 
substituting  0  =  p  into  (2.43)  and  localizing  the  result, 
noting  that  the  volume  integral  of  spatial  mass  density 
over  the  current  configuration  of  the  body  yields  its 
total  mass  which  is  constant  in  time. 

Let  nada  and  N Ad  A  denote  area  elements  on  the 
boundaries  of  m  and  SDT,  related  by  the  usual  Nanson’s 
formula  (Clayton  2014d)  at  a  coincident  point  x  = 
x(X,D,t): 


NAdA  =  jFaAnada.  (2.57) 

Let  dPa  =  ta da  =  tfidA  denote  a  component  of  a  dif¬ 
ferential  mechanical  force  vector,  where  traction  com¬ 
ponents  are  defined  according  to 

ta  =  ctabnb,  tg  =  gab PA  N A .  (2.58) 

The  first  Piola-Kirchhoff  stress  and  Cauchy  stress 
oab  are  related  by  the  deformation  gradient  and  Jaco¬ 
bian  determinant  as  follows  from  (2.58)  and  also  obey 


Po^ 


Ba  +  gabPl 


b\\A * 


(2.61) 


Let  U  be  the  internal  energy  density  measured  per 
unit  reference  volume.  Restricting  the  current  presen¬ 
tation  to  adiabatic  conditions,  i.e.,  no  heat  flux  or  heat 
sources  applied  to  the  body,  the  global  balance  of 
energy  is  stated  in  integral  form  as 

—  f  (—v  v  +  u\dV  =  (i)  *o  • 

DtJm\  2  /  Ism 

+  (£  z  •  DdA  +  f  B  v dV,  (2.62) 

Jdvn  Jm 

where  z  is  a  conjugate  traction-like  force  to  the  time 
derivative  of  internal  state  defined  by 


D(X,t) 


dDA(x,t)  d 
dt  Hd*' 


z  =  zaSDa  =  ZbaNb8Da, 


(2.63) 


with  components  of  a  corresponding  stress-like  ten- 
sor.  The  left  side  of  (2.62)  accounts  for  the  rate  of 
change  of  kinetic  plus  internal  energy,  the  right  side 
for  work  done  by  traction  on  the  boundary  and  by  the 
distributed  body  force.  Use  of  (2.14),  (2.58),  (2.61), 
and  (2.63)  gives  a  global  energy  balance  that  can  sub¬ 
sequently  be  localized  to  material  point  X  as 


u  =  PAv“a  +  {ZaDb)  ||A. 


(2.64) 


The  first  term  on  the  right  is  the  stress  power,  the  second 
the  rate  of  working  from  changes  in  internal  state.  It  is 
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again  noted  that  when  D  is  omitted  from  all  equations, 
the  resulting  balances  of  mass,  momentum,  and  energy 
reduce  to  their  counterparts  from  nonlinear  continuum 
mechanics  (Clayton  2011)  in  Riemannian  geometry 
(and  more  specifically,  Euclidean  space). 


2.2.2  Constitutive  assumptions  and  identities 

Internal  energy  density  per  unit  reference  volume  on 
TO  is  of  the  following  general  form: 


U  =  U(F,  n,  D,  VZ>,  G)  =  U(FaA,  T) ,  Da ,  Dfe,  Gab)- 

(2.65) 


Entropy  density  per  unit  reference  volume  is  the  scalar 
field  rj(X,  D,t).  As  discussed  in  Clayton  (2016a), 
Clayton  (2017a),  physical  reasoning  behind  this  form 
follows  from  generalized  continuum  theories  of  mate¬ 
rials  with  microstructure  including  phase  field  models 
(Capriz  1989;  Clayton  and  Knap  2011a,  2015a;  Levitas 
2014).  The  internal  state  vector  D  is  considered  to  be 
a  vector- valued  set  of  order  parameter(s).  Thermody¬ 
namic  forces  follow  from  the  material  time  derivative 
and  chain  rule  applied  to  (2.65): 


du  .  du . 


+ 


dU  D 
dDA  Dt 


3  T) 

( d\b )  + 


dU  D A 


3  Da 
3  U 


3  GAB 


Gab 


PaK+Ov  +  QAD 


r\A 


+zba2.(dab)  +  sabgab. 


(2.66) 


The  temperature  field  is  0(X,  D,  t).  Spatial  coordinate 
invariance  requires  that  strain  dependence  of  internal 
energy  density  be  of  the  form 


U  =  U[C(F,  g),  rj,  D,  VZ>,  G] 

=  U(CAB,r1,DA,DA,GAB).  (2.67) 

It  follows  that  the  first  Piola-Kirchhoff  stress  PA  and 
Cauchy  stress  oab  obey  a  local  angular  momentum  bal¬ 
ance  that  agrees  with  (2.59): 

pA  r.  Fb  d_£_ 

oCab 

oab  =  jgacPcAFbA  =  2  jFaAFbB22L  =  aba.  (2.68) 

oCab 

Invoking  a  variational  principle  in  prior  work  (Clay¬ 
ton  2017a),  a  local  equilibrium  equation  was  derived  for 


micro-momentum,  containing  terms  in  Q  and  Z.  In  the 
dynamic  regime,  a  rate  equation  for  internal  state  vec¬ 
tor  components  DA  is  posited  by  setting  the  residual  of 
that  equilibrium  equation  proportional  to  the  negative 
rate  of  internal  state: 

DK  =  -L kc[Qc  -  3 AZA  -  ZbHab  +  ZbHac 

—  dBZAdADB  —  ZB(dcNB 

-3  cKadDd  +8aC°dDbb) 

-  pAdBdc<padADB 

+  ( SAB  +  UGAB)dcGAB].  (2.69) 

Here,  LKC  is  a  positive  definite  matrix  of  constants 
(depending  on  the  material)  that  control  the  time  scale 
for  the  rate  of  change  of  internal  state.  Equation  (2.69) 
states  that  the  order  parameters  evolve  in  time  such  that 
at  equilibrium,  the  term  in  square  braces  vanishes  in 
accordance  with  the  static  director  momentum  equation 
derived  previously  (Clayton  2017a).  Kinetic  equation 
(2.69)  can  be  expressed  in  condensed  form  as 


bK 


|KC 


"  du  /  du  \ 
Jbc  ~  Va  vavADc ) 


+  •••  , 

(2.70) 


where  here  V  is  the  horizontal  covariant  derivative  and 
remaining  terms  have  been  truncated  only  for  presenta¬ 
tion  purposes.  This  equation  is  reminiscent  of  the  time 
dependent  Ginzburg-Landau  or  Allen-Cahn  equations 
(Allen  and  Cahn  1979;  Levitas  2014). 

Several  choices  depending  on  the  class  of  material 
complete  the  model.  A  metric  tensor  G  is  introduced, 
from  which  connection  coefficients  are  obtained  from 
relations  in  Sect.  2.1.1.  Horizontal  and  vertical  con¬ 
nection  coefficients  in  (2.8)  and  (2.9)  must  be  chosen, 
typically  those  of  the  Chern-Rund  connection.  Simi¬ 
lar  features  are  assigned  for  the  current  configuration, 
as  described  by  equations  of  Sect.  2.1.2.  The  internal 
energy  density  function  U  in  (2.65)  must  be  assigned. 
Finally,  constitutive  equations  for  inelastic  components 
of  the  deformation  gradient,  FD  of  Sect.  2.1.5,  are  usu¬ 
ally  needed,  as  are  kinetic  coefficients  LKC  for  general 
dynamic  problems. 


2.2.3  Multiplicative  thermodynamics 

Extending  prior  treatments  of  Clayton  (20 1 6a,  2017  a,  b) 
to  account  for  entropy,  and  those  addressed  in  phase 
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field  theory  (Clayton  and  Knap  2011a,  2015a),  the 
internal  energy  U  of  (2.65)  is  split  into  a  sum  of  ther¬ 
moelastic  energy  density  W  and  structure/intemal  state 
dependent  energy  density  /: 

U{F%,  r,,  Da,  Da,  Gab )  =  W[(FE)“,  r1,  DA,  gafi] 

+f(DA,DA,GAB). 

(2.71) 

Notice  that  general  dependence  on  F  is  replaced  by 
dependence  on  its  thermoelastic  part  FE .  A  more  gen¬ 
eral  treatment  would  allow  for  couplings  among  elas¬ 
tic  deformation,  entropy,  and  internal  state.  For  terms 
entering  the  multiplicative  kinematics,  the  conditions 
(2.44)  and  (2.45)  are  imposed: 

(FEra[F,  Fd(D )] 

=  FaA[X,  D(X,  t),  tWD~l)aW{X,  t )].  (2.72) 

Applying  (2.72)  and  invoking  (2.49)  for  intermediate 
space  metric  components  gap,  thermodynamic  forces 
of  (2.66),  at  fixed  X,  obey 


du 

dF° 


dW 

JfI 


e  = 


Qa  = 


yB  _ 


dU  _ 
dr]  dr] 
dU  dW 


dDA 

dU 

jDfB 


(FD~X 


9/ 


dDA 

9/ 

dDAB 


dW 


b,vE^FDTb 


+  ^-pa(ptra 


dDA 


(2.73) 


of  a  quantity  upstream  and  downstream  from  the  shock. 
The  jump  of  a  quantity  across  the  shock  plane  is  then 

[•]  =  (•)“-  (•)+•  (3.1) 

Denote  by  w  =  v  —  D  the  velocity  of  the  material  rela¬ 
tive  to  the  shock  front,  with  v  the  scalar  particle  veloc¬ 
ity,  i.e.,  the  component  of  velocity  vector  v  parallel  to 
the  direction  of  shock  propagation.  For  convenience, 
let  the  shock  direction  be  parallel  to  X1,  and  further 
assume  that  uniaxial  strain  conditions  hold: 

F  =  1  +  (du/dX)g1  (g)  G1 .  (3.2) 

Basis  vectors  g1  and  G1  point  in  the  direction  of  motion 
X  =  X1,  and  u(X,  D,  t)  is  the  displacement  compo¬ 
nent  in  this  direction.  The  equal  first  Piola-Kirchhoff 
and  Cauchy  stress  components  normal  to  the  front  are 
equal  to  the  negative  of  the  so-called  shock  pressure, 
i.e.,  =  cr l  =  —P.  The  internal  energy  per  unit  mass 
is  e  =  U/po. 

The  form  of  Reynolds’  transport  theorem  derived 
in  (2.43)  can  be  used  in  conjunction  with  assumptions 
and  definitions  for  time  derivatives  in  Sect.  2.1.4  to 
obtain  a  general  transport  theorem  accounting  for  flux 
of  a  quantity  (j)(x,  D,t)  across  the  shock  plane.  The 
derivations  of  Casey  (Casey  2011)  (in  Euclidean  space) 
apply  in  full  when  the  above  definitions  and  caveats  for 
Finsler  space  hold.  Specifically,  let  subscripts  1  and  2 
denote  subregions  of  the  body  manifold  partitioned 
by  a  shock  plane  X  (, t ).  Then  the  transport  theorem  can 
be  expressed  as  (Casey  2011) 


Furthermore,  dU/dGAB  dU/dGAB  and  Gab(X) 
=  0,  so  SAB  ->  0  in  (2.66)  and  (2.69)  without  energetic 
consequence.  Spatial  invariance  analogous  to  (2.67) 
follows  from  forcing  W  to  depend  on  C£  rather  than 
Fe  ,  where 

0 CEU  =  (FE)aJab(FE)bp  =  (FD-l)aCAB(FD-l)p , 
Cab  =  FaAgabFbB.  (2.74) 


0  =  —  f  ct>dv  =  0i  +  02  +  f  Uwjda.  (3.3) 
Dt Jm  Js 

The  first  two  terms  following  the  second  equality 
account  for  rates  of  change  of  0  in  each  subregion. 
First  taking  0  =  p  in  (3.3),  the  global  form  of  mass 
conservation  law  (2.56)  requires  that  0  vanish  for  each 
subregion  as  well  as  for  the  whole  body,  leaving 


3  Theory:  moving  surfaces  of  discontinuity 

Consider  a  planar  shock  moving  with  a  normal  compo¬ 
nent  of  natural  velocity  D  in  the  material  manifold  3DT, 
across  which  velocity,  stress,  and  deformation  gradient 
may  be  discontinuous.  Let  (-)+  and  (-)_  denote  values 


[pw]  =  0.  (3.4) 

Next  taking  0  =  pv  as  the  spatial  linear  momentum, 
conservation  law  (2.60)  requires 

{P  +  pvw]  =  0.  (3.5) 
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Let  f  denote  the  Piola  transform  of  z  in  (2.62): 


t;Ada  =  ZAnada  =  j  F^ZAnada  =  ZANBdA  =  z^dA. 

(3.6) 

Then  the  conservation  law  for  energy  is  obtained  by 
setting  0  =  p(e  +  ^v2)  in  (3.3)  and  using  (2.62)  and 
(3.6): 

lp\N(e  +  v2/2)  +  Pv  -  ?  •  I>]  =  0.  (3.7) 

Equations  (3.4)  and  (3.5)  are  identical  in  form  to 
the  classical  Rankine-Hugoniot  conditions,  while  (3.7) 
reduces  to  its  classical  counterpart  when  f  •  D  =  0. 

The  equations  simplify  further  when  the  material  is 
undeformed,  unstressed,  and  at  rest  ahead  of  the  shock, 
with  null  internal  energy  (f/+  =  0)  taken  as  the  datum. 
In  that  case,  algebraic  manipulations  yield 

du 

V  =  -D— ,  P  =  povD, 
dX 

U  =  ip0v2  -  2[?  •  DJ.  (3.8) 

Quantities  in  each  equality  correspond  to  the  down¬ 
stream  (shocked)  state  besides  the  jump  term  on  the 
right  side  of  final  equation  in  (3.8).  The  latter  term, 
written  in  general  vector  form  for  the  moment,  is  sub¬ 
jected  to  further  assumptions  later  upon  specification 
of  a  particular  constitutive  model. 


4  Application:  shock  loading  of  boron  carbide 

The  dynamic  theory  of  Finsler-geometric  continuum 
mechanics  is  now  applied  to  shock  loading  of  ceramic 
single  crystals.  In  particular,  the  problem  considered  in 
Sect.  4  is  focused  on  anisotropic  boron  carbide  (B4C) 
subjected  to  planar  impact  loading  along  its  c-axis 
([0001]).  The  pseudo-Finsler  metric  quantifies  local 
densification  (i.e.,  local  volume  decrease)  that  accom¬ 
panies  the  stress-induced  phase  transformation  from  a 
trigonal  crystal  to  glassy  phase  (Yan  et  al.  2009;  Clay¬ 
ton  2012a,  2013a,  2014a;  Clayton  and  Tonge  2015; 
Taylor  et  al.  2012;  Taylor  2015;  An  and  Goddard 
2015a).  Inelasticity  arises  from  two  contributions:  the 
aforementioned  volume  change  and  shear  strain  within 
amorphous  bands  (Chen  et  al.  2003;  An  and  Goddard 
2015a).  The  boundary  value  problem  involves  steady 
shock  wave  propagation,  under  conditions  of  uniaxial 
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strain  compression,  of  a  nonlinear  elastic  domain  of  the 
ceramic.  The  presentation  that  follows  proceeds  from 
a  general  constitutive  description  to  a  focused  solution 
of  the  specific  problem. 


4.1  Geometry  and  deformation 


The  material  body  is  a  domain  of  effectively  infinite 
dimensions.  Uniaxial  strain  conditions  are  imposed, 
with  compression  in  the  X3 -direction.  The  internal 
state  vector  is  prescribed  as  {DA}  {0,  0,  D3},  later 
linked  physically  to  structure  collapse  of  boron  car¬ 
bide.  By  construction,  field  variables  may  vary  only 
with  Z  =  X3  and  D  =  D3.  Components  Dl  and  D 2 
are  superfluous,  while  X1  and  X2  are  important  in  the 
context  of  lateral,  i.e.,  transverse,  stresses,  for  example. 
A  Cartesian  coordinate  system  is  used  for  { XA },  and 
thus  the  metric  tensor  G  contains  no  explicit  depen¬ 
dence  on  X. 

Definitions  and  identities  of  Sect.  2. 1  are  given  more 
specific  forms  to  be  consistent  with  these  protocols.  An 
isotropic  metric  tensor  is  prescribed,  similar  to  that  in 
prior  work  (Clayton  2016b): 

{X,  Y,  Z}  =  {X\X2,X3},  D  =  D3  =  D(Z,t); 


G(D)  =  G(D)  =  B(D)  1  = 


(4.1) 

~B(D )  0  0 

0  B(D)  0 
0  0  B(D)_ 

Gab  =  BSab ;  G  =  det  G  =  B3(D);  (4.2) 

Yabc  =  2(^aGbc  +  3 bGac  ~  3 cGab)  =  0, 

Ga  =  \ybCDBDc  =0,  N§  =  dBGA  =  0 

=►  3a(-)  =  3aO;  (4.3) 


C131  =  C311  =  C333  =  C322  =  C232  =  B'/2, 

Ci  13  =  C223  =  ~B'/2.  (4.4) 

The  prime  notation  obeys  B'(D )  =  dB(D)/dD ,  and 
all  other  covariant  components  Cabc  of  Cartan’s  tensor 
are  zero.  The  Cartesian  metric  1  is  scaled  isotropically 
by  scalar  function  B ,  a  conformal  transformation  of 
Weyl  type  (Weyl  1952;  Clayton  2016a,  2017a).  The 
pseudo-Finsler  manifold  DJI  is  locally  Minkowskian 
Minguzzi  (2014). 

Components  of  coefficients  of  the  Chern-Rund  con¬ 
nection  are  used,  also  following  prior  work  (Clayton 
2016b,  2017a).  From  this  choice  and  from  the  vanish¬ 
ing  nonlinear  connection  coefficients  in  (4.3),  it  follows 
that 


Finsler-geometric  continuum  dynamics  and  shock  compression 


65 


11 A  A  r~f  r4 

nBC  ~  KBC  ~  1  BC 

=  \Gad  (ScG  bd  +  SbGcd  —  $dGbc)  =  0; 
Vbc  =  Y'bc  =  0-  (4-5) 

For  the  spatial  configuration,  i.e.,  the  deformed 
material  manifold  m,  coordinates  and  metric  tensor  are 
of  fully  analogous  forms  to  those  used  for  TO.  Invoking 
lower  case  font  for  the  spatial  frame, 


Motions,  deformations,  and  director  gradients 
defined  in  Sect.  2.1.3  take  the  following  forms  under 
unaxial  strain  conditions.  Let  z  =  cp(Z,  D,t)  denote 
deformation  in  the  Z  direction  with  e(Z,  D,t)  a  cor¬ 
responding  displacement  gradient  measure,  negative  in 
compression,  such  that 

x  =  X,  y  =  Y,  z  =  (p\ 

d  =  d(Z,D,t );  D  =  D(Z,t );  (4.10) 


F(X,  D,  t) 


dx(X)/dx  dx(X)/dY  dx(X)/dz 

dy(Y)/dX  dy(Y)/dY  dy{Y)/dZ 

d(p(z ,  d ,  t)/dx  d<p(z ,  d ,  t)/dY  dy(z,  d ,  t)/dz 


T 

0 

0 

"1 

0 

0 

0 

1 

0 

= 

0 

1 

0 

0 

0 

d<p(Z,  D,t)/8Z_ 

_0 

0 

1  +e(Z,  D,t)_ 

(4.11) 


{x,  y,  z}  =  {xl,x2,  x3},  d  =  d3  =  d(z,  t ); 
g(d)  =  g(d)  =  bid)  1, 

g(d)  =  b\d );  (4.6) 

K<ji>c  =  9  (.^agbc  +  dbgac  dcgab)  =  0, 

=  0, 

nab  =  dbga  =  0  =►  «fl(.)  =  3fl(0;  (4.7) 

C131  =  C311  =  C333  =  C322  =  C232  =  y  /i, 

C113  =  C223  =  (4.8) 

where  Z/(d)  =  db(d)/dd.  The  spatial  geometry  is  also 
locally  Minkowskian.  Using  Chern-Rund  connection 
coefficients  with  vanishing  nonlinear  coefficients  from 
(4.7), 

it  cl  j^ci  r^a 

nbc  —  Kbc  —  1  be 

=  \gad(&cgbd  +  hged  -  &dgbc)  =  0; 

Vbc  =  Ybc  =  0-  (4.9) 

Summarizing,  nonlinear  connection  coefficients 
(Ng  and  n^)  vanish  identically  in  both  configurations, 
as  do  Chern-Rund  coefficients  (r£c  and  rgc).  Cartan’s 
coefficients  C^c  and  C£c  may  be  nonzero.  The  hori¬ 
zontal  covariant  derivatives  of  the  metric  tensors  van¬ 
ish  so  (VG)|A  =  0  and  C \/g)\a  =  0-  Thus,  use  of 
Rund’s  version  of  Stokes’  theorem  of  (2.14),  (2.15), 
and  (2.25)  to  derive  balance  and  transport  equations  is 
mathematically  valid.  Vertical  connection  coefficients 
(V BC^bc^BC^  anc*  ^c)  all  vanish  by  definition,  which 
greatly  simplifies  calculations. 


J  =  F\FlFl4ifG  =  Fi(g/G )J/2 

=  (1  +e)(^/B)3/2;  (4.12) 

D\ 3  =  d3D  -  yV|  +  K^D  =  dD/dZ  =  D' ,  (4.13) 

where  the  prime  notation  in  (4.13)  denotes  a  partial 
derivative  with  respect  to  the  axial  coordinate. 

The  non- vanishing  component  of  the  director  vector, 
i.e.,  the  internal  state  variable  D ,  is  physically  related 
to  transformation  to  a  glassy  phase  and  shear  failure  of 
boron  carbide.  An  overall  loss  of  shear  strength  and 
volume  collapse  are  simultaneously  associated  with 
nonzero  values  of  D.  A  regularization  constant  l  with 
dimensions  of  length  and  a  normalized  order  parameter 
§  G  [0,  1]  are  now  introduced: 

§  =  D/l  £'  =  D'/l.  (4.14) 

The  scalar  /  is  identified  as  the  value  of  state  variable 
D  at  which  the  material  ruptures  and  undergoes  com¬ 
plete  densification,  and  it  is  assumed  to  be  a  material 
property.  Letting  k  denote  a  constant  depending  on  the 
material,  a  more  specific  form  of  the  pseudo-Finsler 
material  metric  in  (4.2)  is  now  prescribed: 

G(D)  =  B(D)  1  =  exp[(£/3)(£>/02]l; 

G(0  =  B 3(£)  =  exp(^2) 

=>■  3B'/(2B)  —  kD/l  =  £§.  (4.15) 

The  third/axial  component  of  the  trace  of  Cartan’s  ten¬ 
sor  of  (4.4)  will  be  used  later: 

C3A  =  Cj  C3AB  =  U  C311+U  C322  +  Cj  C333 

=  3B'/(2B)  =  kD/l  =  k%.  (4.16) 
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This  form  of  conformal  transformation  is  similar,  but 
not  identical,  to  that  used  in  Clayton  (2017a).  A  ben¬ 
efit  of  the  current  prescription  is  vanishing  of  (4.16) 
at  §  =0,  which  eliminates  certain  driving  force  con¬ 
tributions  for  transformation  in  the  fully  elastic  state. 
The  length  of  a  referential  line  element  in  (2.10)  and 
the  corresponding  volume  form  in  (2.11)  become 

|dX|2  =  exp[(£/3)(£>//)2](dX  •  dX 
+  dY  ■  dY  +  dZ  ■  dZ), 

dtt  =  exp[(k/2)(D/lf]dX  A  dY  A  dZ.  (4.17) 

For  §  >  0,  expansion  occurs  when  k  >  0  and  contrac¬ 
tion  when  k  <  0.  Volume  collapse  in  boron  carbide  is 
associated  with  the  latter  condition. 

As  in  prior  work  (Clayton  2016b,  2017a,  b),  the  spa¬ 
tial  and  referential  state  vector  entries  are  set  to  coincide 
in  the  sense  (4.10): 

d(Z,  D,  t)  =  D{Z ,  t)  =  /§(Z,  t).  (4.18) 

This  can  be  interpreted  as  a  scalar  push-forward  relation 
for  field  D  (Clayton  2017b).  An  analogous  Weyl-type 
form  of  the  spatial  metric  of  (4.6)  is 

g(d )  =  b(d)  1  =  txp[(k /3)(D  / Z)2]l; 

*(0  =  b\$)  =  exp  (£|2)  =►  3b' /(2b)  =  kD/l  =  U-. 

(4.19) 

With  these  choices  of  metrics,  (4.12)  reduces  to 


J(Z,  D,  t )  =  F33(Z,  £>,0  =  1+  C(Z,  D,  t).  (4.20) 

As  remarked  already,  inelasticity  in  shock  com¬ 
pressed  boron  carbide  consists  of  contributions  of  den- 
sification  as  well  as  shearing  in  amorphous  bands.  For 
compressive  loading  normal  to  the  c-axis,  densifica- 
tion  is  presumed  spherical  as  in  Clayton  (2016b).  Also, 
for  loading  according  to  this  protocol,  inelastic  shear¬ 
ing  takes  place  on  planes  and  in  directions  analogous 
to  those  for  pyramidal  slip  of  the  type  (1 101)  {01 1 1} 
in  hexagonal  crystals,  following  examination  of  recov¬ 
ered  fragments  (Chen  et  al.  2003;  Yan  et  al.  2009). 
Shear  localization  and  amorphization  for  this  kind  of 
deformation  system  were  studied  via  atomic  simula¬ 
tions  of  simple  shearing  of  boron  carbide  single  crys¬ 
tals  in  An  and  Goddard  (20 15a),  wherein  no  twinning  or 
dislocation  glide  were  reported.  These  authors  further 
considered  shearing  of  a  different  crystal  orientation 
(1010){0001|  resulting  in  twinning  and  amorphization 
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on  the  basal  plane.  The  present  study  does  not  address 
basal  plane  localization  or  twinning  since  elastic  driv¬ 
ing  forces  (i.e.,  resolved  shear  stress)  on  (0001)  vanish 
and  since  such  modes  are  inconsistent  with  symmetry 
for  c-axis  compression  of  the  crystal. 

The  multiplicative  description  of  (2.44)  is  used,  with 
the  inelastic  portion  containing  micro  structure  (i.e.,  D) 
dependent  terms.  The  total  deformation  gradient  in 
Cartesian  coordinates  is 


F  = 


1  0 
0  1 
0  0 

r(FEyx 

o 
o 

r{FD)l 

0 

0 


0 

0 

1+6 


=  fe  fl 


0 

(FE)yY 

0 

0 

(■ FD)Yy 
0 


0 

0 

(■ FE)ZZ 

0 

0 

D\Z 


(F  ) 


z  J 


(4.21) 


The  inelastic  term  is  diagonal  due  to  symmetry  of 
the  crystal  structure  and  the  present  loading  mode.  It 
accounts  for  the  inelastic  volume  change  and  inelas¬ 
tic  shearing  yo  on  the  six  pyramidal  systems  of  type 

<1101>{01ll}: 


(4.22) 


Here,  the  material  constant  x  =  ^[exp(k/2)  —  1] 
accounts  for  isotropic  volume  collapse  with  k  <  0. 
Unit  vectors  sa  and  ma  are  the  orthogonal  direction  and 
plane  for  inelastic  pyramidal  shear  of  maximum  mag¬ 
nitude  yo,  another  material  constant.  The  terms  con¬ 
taining  yo  are  the  first  three  in  the  series  approximation 
of  the  exponential  function  corresponding  to  an  exact 
representation  of  this  symmetric  inelastic  deformation 
mode  (Clayton  2014b,  e).  Finally,  *(§)  :  [0,  1]  -> 
[0,  1]  is  an  interpolation  function  with  vanishing  end¬ 
point  derivatives,  specifically  as  used  in  phase  field  rep¬ 
resentations  of  structural  changes  (Levitas  et  al.  2009; 
Clayton  and  Knap  2011a): 
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t(£)  =  3£2-2£3,  t'($)  =  6|(1-|).  (4.23) 

The  form  of  (4.22)  scales  contributions  of  phase  trans¬ 
formation  and  inelastic  shearing  by  a  phase  field- type  of 
interpolation  function  (4.23).  The  first  term  in  square 
braces  is  assumed  isotropic  or  spherical,  accounting 
for  compaction  commensurate  with  the  crystal-to-glass 
transition.  Depending  on  the  loading  protocol  and 
deformation  system,  the  phase  change  could  demon¬ 
strate  some  anisotropy  (e.g.,  preferential  structure  col¬ 
lapse  along  the  c-axis  (Yan  et  al.  2009;  Clayton  2012a)), 
but  the  present  work,  as  in  Clayton  (2016b),  simplifies 
this  part  of  transformation  strain  as  isotropic.  Terms 
involving  the  dyadic  product  of  shearing  direction  and 
plane  encompass  all  deviatoric  mechanisms  associated 
with  amorphization,  shear  localization,  and  eventual 
mode  II  failure  on  pyramidal  planes.  The  assumption 
that  evolution  of  volumetric  and  deviatoric  inelastic 
deformations  proceed  in  lock-step,  via  the  same  func¬ 
tion  l  ,  benefits  from  simplicity  and  is  demonstrated  later 
to  depict  physically  reasonable  solutions  to  the  present 
highly  symmetric,  one-dimensional  shock  problem  for 
compressive  strains  up  to  on  the  order  of  20%.  Gen¬ 
eralizations,  albeit  at  the  possible  expense  of  further 
assumptions  or  calibration,  are  anticipated  to  be  neces¬ 
sary  in  the  future  to  address  more  sophisticated  prob¬ 
lems. 

Noting  that  slip/shear  is  isochoric,  the  inelastic  vol¬ 
ume  change  is  approximately 


strain  tensor  more  accurately  represents  the  nonlinear 
elastic  response  of  strong  solids  with  a  large  ratio  of 
effective  shear  to  bulk  modulus  such  as  quartz  (Clay¬ 
ton  2014b,  2015c)  and  boron  carbide  (Clayton  and 
Tonge  2015)  than  the  Green  elastic  strain  (Clayton  and 
McDowell  2003;  Clayton  2011)  or  Eulerian  material 
strain  (Clayton  2013b;  Lloyd  et  al.  2014b).  The  loga¬ 
rithm  of  the  elastic  stretch  corresponding  to  the  first  of 
(2.74)  is 

eE  =  In  UE  =  ln[(C£)1/2]  =  1  In  CE, 

(eE)ap  =  ^(\nCE)ap.  (4.25) 

For  the  present  class  of  problems  with  (4. 1 1)  and  (4.21) 
now  invoked  and  gap  =  8ap,  the  first  of  (2.74)  results 
in  the  following  three  possibly  nonzero  components: 

(cE) !  =  [( feyx ]2,  ( cE)\  =  [c FEyY]2 , 

0 CE)l  =  [(FE)zzf.  (4.26) 

Computation  of  the  thermoelastic  logarithmic  strain 
(omitting  the  redundant  numerical  superscript)  is  then 
trivial  since  CE  is  diagonal: 

ef  =  (. eE)l  =  \n[(FE)xx], 
el  =  (, eE)\  =  ln[(FE)yY], 

eE  =  (eE)l  =  ln[(FE)zz].  (4.27) 

Thermoelastic  volume  change  is  related  to  the  trace  of 
the  logarithmic  strain  via  In  JE  =  ef  +  ef  +  ef . 


JD  =  det  Fd  ~  (1  +  ^x)3  ~  1  +  3^x  =  £exp(&/2). 

(4.24) 

Crystal  plasticity-type  representations  of  inelastic 
shearing  modes  distinct  from  slip  have  been  used  else¬ 
where  for  cleavage  fracture  in  rocks  (Clayton  2010) 
and  other  brittle  solids  (Aslan  et  al.  2011). 

Parameter  yo  is  regarded  as  a  constant  for  simplic¬ 
ity,  here  describing  shearing  on  pyramidal  planes.  For 
other  deformation  systems  or  other  shear  mechanisms 
in  boron  carbide,  yo  would  presumably  vary,  and  it 
would  also  presumably  vary  with  material  composi¬ 
tion.  In  (4.22),  function  i(t-)  scales  the  magnitude  of 
shear  deformation,  where  evolution  of  §  is  dictated  by  a 
kinetic  process  or  incremental  equilibrium  conditions. 

The  material  logarithmic  thermoelastic  strain  tensor 
eE  is  used  in  the  constitutive  model  of  boron  carbide, 
as  in  Clayton  and  Tonge  (2015),  Clayton  (2017b).  This 


4.2  Thermomechanics 


The  internal  energy  per  unit  reference  volume  U  in 
(2.65)  is  additively  split  into  a  thermoelastic  strain 
energy  per  unit  reference  volume  W  and  a  structure 
or  phase  dependent  energy  per  unit  reference  volume 
/,  following  (2.71).  The  following  sum  is  used: 


U{F%,  r,,  Da ,  DAB,  GAB)  =  W(eEp,  rj)+f(DA,  DAB); 

(4.28) 

W  =  +  \c^e%eEyie% 

+  0qt,[1  -  rfe%  +  t]/(2co)];  (4.29) 


r 


r 


/  =  7|vz>l  +/3|Z)|' 


(4.30) 


The  thermoelastic  strain  energy  function  W  is  essen¬ 
tially  identical  to  one  used  in  Clayton  and  Tonge  (2015), 
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accounting  for  anisotropic  linear  and  nonlinear  elas¬ 
tic  effects,  thermoelastic  coupling,  and  specific  heat. 
Entropy  r)  —  0  and  temperature  0  =  0q  for  the  datum 
of  the  ambient  reference  state.  Isentropic  second-  and 
third-order  elastic  constants  are  Ca^y8  and  Ca^y8e(^ . 
These  are  taken  herein  to  be  independent  the  phase  of 
the  solid  since  inelastic  deformation  rather  than  mod¬ 
uli  degradation  as  in  Clayton  (2016b)  now  accounts 
for  shear  softening.  The  specific  heat  at  constant  elas¬ 
tic  strain  per  unit  reference  volume  is  the  constant 
co.  The  anisotropic  and  symmetric  Griineisen  tensor 
most  generally  consists  of  six  constants  T0  .  Since 
Gab(D)  =  8abB(D ),  explicit  dependence  of  func¬ 
tions  on  the  right  side  of  (4.28)  on  the  metric  tensor 
would  be  redundant. 

The  microstructure  dependent  function  /  consists  of 
three  parts:  a  form  quadratic  in  \D\,  representative  of 
fracture  or  rupture;  a  form  quadratic  in  |  VZ>  | ,  account¬ 
ing  for  energy  of  phase  boundaries,  shear  bands,  and/or 
crack  surfaces;  and  a  double  well  potential  of  order  four 
in  \D\.  The  surface  energy  per  unit  reference  area  is 
the  intrinsic  material  constant  T.  This  energy  can  be 
associated  with  the  formation  of  mode  II  shear  zones 
commensurate  with  cavitation  and  rupture  upon  com¬ 
plete  amorphization  (An  and  Goddard  2015a).  The  first 
term  on  the  right  side  of  (4.30)  is  the  standard  quadratic 
form  for  gradients  of  order  parameter(s)  in  models  of 
fracture  (Clayton  and  Knap  2014,  2015b)  and  phase 
transformations,  e.g.,  Levitas  (2014).  The  second  term 
on  the  right  side  of  (4.30)  is  the  usual  prescription  for 
bulk  fracture  energy  in  phase  field  theory  (Clayton  and 
Knap  2014,  2015a).  Finally,  the  double  well  poten¬ 
tial  with  barrier  strength  quantified  by  A,  a  constant 
of  dimensions  of  energy  per  unit  reference  volume, 
is  the  conventional  form  used  in  phase  field  models 
of  phase  transformations  and  twinning  (Clayton  and 
Knap  2011a;  Levitas  2014).  In  the  present  work,  the 
phase  transformation  corresponding  to  the  double  well 
potential  is  that  from  crystal  to  amorphous  solid  in  B4C. 
In  summary,  /  is  a  hybrid  potential  combining  standard 
forms  from  the  literature  to  account  for  energies  asso¬ 
ciated  with  related  processes  of  phase  transformation, 
shear  localization,  and  shear  failure  in  shock  loaded 
boron  carbide  crystals. 

The  generic  energy  densities  W  and  /  are  applica¬ 
ble  to  any  loading  direction  and  any  Cartesian  coor¬ 
dinate  system,  with  full  anisotropy  and  a  vector¬ 
valued  internal  state  variable,  respectively.  These  gen¬ 
eral  forms  are  simplified  further  for  the  present  appli- 
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cation  in  shock  physics  with  an  effectively  scalar  order 
parameter.  Regarding  thermoelastic  strain  energy,  the 
Griineisen  scalar  is  used  for  thermoelastic  coupling, 
and  an  isotropic  contribution  from  third-order  elasticity 
is  used  to  account  for  increases  in  tangent  bulk  modulus 
with  increasing  mass  density  under  shock  compression. 
Letting  Bo  and  B'0  denote  the  isentropic  bulk  modulus 
and  its  pressure  derivative  at  the  reference  state,  and 
letting  r^j  — >►  7~b  the  strain  energy  function  in 
(4.29)  becomes 

w  =  l-C uefej  +  Iflo(2  -  B’0)( In  JE )3 

+  d0r][l-  r0lnJE  +  r1/(2c0)].  (4.31) 

Indices  I,  J  =  1,2,  ...6  denote  Voigt  notation,  and 
C  ij  are  the  second-order  elastic  constants  to  be  intro¬ 
duced  for  boron  carbide  in  Sect.  4.3.  Now  considering 
(4.13)  and  (4.14)  whereby  D(Z,  t)  =  /§(Z,  t)  is  the 
only  relevant  component  of  D,  the  phase  transition  and 
fracture  energy  function  in  (4.30)  becomes 

/  =  Tl[$2/l2  +  (£')2]  +  Af(  1  -  £)2.  (4.32) 


More  specific  forms  of  thermodynamic  forces  for 
Finsler  continuum  mechanics  at  a  fixed  material  point 
X  are  then  obtained: 
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(4.33) 

(4.34) 


(4.35) 

(4.36) 


For  the  present  problem  of  shock  compression,  now 
consider  the  Cartesian  coordinate  system  and  uniax¬ 
ial  strain  conditions  of  Sect.  4.1.  Boron  carbide  single 
crystals  have  rhombohedral,  i.e.,  trigonal  symmetry. 
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The  current  application  assigns  the  XT-plane  as  the 
basal  plane,  i.e.,  (0001)  in  hexagonal  Miller-Bravais 
indices,  with  Z  the  direction  of  loading  along  [0001]. 
For  this  loading  mode,  pertinent  isentropic  second- 
order  elastic  constants  in  Voigt  notation  are  C 1 1  =  C22 , 
C33,  C12,  and  C13  =  C23.  Other  loading  directions 
or  more  general  crystal  orientations  might  require  use 
of  C14  =  —  C24  =  C56  which  are  not  needed  here. 
The  thermoelastic  strain  energy  in  (4.29)  and  (4.31) 
becomes,  for  the  potentially  nonzero  elastic  strain  com¬ 
ponents  in  (4.27), 


dFD(l=) 


=  *'($) 


Xl  +  yo^s01 


a=l 


2/6 


\a=l 


\a=\  / 


(4.41) 


For  shock  loading  in  the  Z  direction,  define  the 
shock  stress,  positive  in  compression,  as 


w  =  W(e f ,  el ,  el,  /?)  =  l[Cn(ef)2  +  Cn(ef)2 

+  C33(ef)2  +  2(Ci2«f«f 

+  Ci3efe^  +  C\3eEeE)\ 

+  .50(2  —  Bo)(ef  +  e2  4”  e3)2 

o 

+  0oriU  ~  r0(e f  +  el  +  e3£)  +  r,/(2c0)l 

(4.37) 

For  the  present  loading  mode  and  symmetry  consider¬ 
ations,  lateral  elastic  strains  are  equal,  i.e., 

(FE)xx  =  (FE)yY,  ef  =  ef.  (4.38) 

Three  possibly  nonzero  stress  components,  two  of 
which  are  equal,  are  obtained  from  (4.33)  and  (4.37): 

Pl  =  Px  =  P2  =  Py  =  (Cll  +  Ci2)ef 

+  Cl3el  +  1b0(2  -  So)(ln  JE)2  ~  Oor0ri, 

P]  =  PzZ  =  ^[20^+0^ 
+^B0(,2-Bl))(lnJE)2-e0r0r]  . 

(4.39) 


The  conjugate  thermodynamic  force  to  D  =  Z£  in 
(4.35)  is  then  given  by 


Q  =  -  7 


2y£  +  2A§(1  —  £)(1  —  2§) 


,  F  x  d(FD)* 

-2P{(FEyx^j^x 


-  p3(fe)zz 


d(Fu) 


D\Z 


(4.40) 


where  from  (4.22), 


P  =  -P33  =  -Pf .  (4.42) 

The  treatment  of  surfaces  of  discontinuity  of  Sect.  3 
applies  here  if  direction  X  is  replaced  with  direction 
Z.  The  thermodynamic  driving  force  in  (3.6)  is,  for 
uniaxial  strain  and  a  single  scalar  component  of  Z>, 

f  =  ft  =  T  =  2  T£'.  (4.43) 

For  the  present  problem,  the  jump  of  scalar  product 
(  D  across  the  shock  front  becomes 

[?  -i>]  =  2r /[$'£].  (4.44) 

Regions  far  ahead  and  far  behind  the  shock  front  are 
assumed  to  be  in  equilibrium  with  respect  to  internal 
state.  Across  the  front  of  presumed  width  /,  the  fol¬ 
lowing  diffuse  interface  approximations  are  imposed 
for  the  jump  in  order  parameter  gradient  and  its  rate, 
whereby  §  increases  from  zero  (i.e.,  its  far  upstream 
value  §+  =  0)  to  its  downstream  value  with 
decreasing  material  coordinate  Z : 

=►  2771^1/D  «  -2(77 Z)((T)2.  (4.45) 

Now  let  — >  §  in  subsequent  equations  since  §  van¬ 
ishes  upstream.  Analogs  of  Rankine-Hugoniot  equa¬ 
tions  in  (3.8)  then  become 

v  =  — De,  P  =  p0vD,  U  =  ip0v2  +  2y£2. 

(4.46) 

Downstream,  the  Ginzburg-Landau  type  equation 
(2.69)  presumably  holds  with  vanishing  left  side,  i.e., 

D  =  lk=  0. 

Forthcoming  equations  pertain  to  the  downstream 
state.  The  shocked  material  in  this  state  has  nonzero 
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particle  velocity  v  and  density  p,  but  it  has  null  accel¬ 
eration.  Stress  does  not  vary  with  time  in  regions  far 
behind  the  shock  so  axial  stress  is  of  functional  form 
P  =  P(Z,  D).  The  only  significant  equation  in  lin¬ 
ear  momentum  balance  (2.61)  describes  the  direction 
of  compressive  loading:  a  =  z  =  3.  Recalling  that 
off-diagonal  components  of  P  vanish,  momentum  con¬ 
servation  in  the  downstream  state  requires  P^  =  0. 
Substituting  from  (2.15)  and  (4.16),  this  results  in  the 
ordinary  differential  equation 


dP{Z ,  D)  dP(Z ,  D)  dD  ^ 


az 

d  P 

=  —  +  P 
dZ 


dD 
kD  d D 
~TdZ 


az 
=  o. 


3B\D)  dD 
- —  P - 

2  b(d)  az 


(4.47) 


The  other  two  macroscopic  linear  momentum  balance 
equations  simply  require  that  transverse  stress  com¬ 
ponents  P*  =  Pj  are  constant  with  respect  to  X 
and  Y,  respectively.  These  requirements  are  consistent 
with  the  assumed  X  and  Y  independence  of  solution 
fields  for  unaxial  loading.  The  only  pertinent  equation 
in  kinetic  law  (2.69)  is,  with  §  =  0  for  the  downstream 
state, 


ar(z,  D) 


az 


+ 


ar(z,  D)  +  3b\d)t 


-p 


d2p{Z ,  D) 
dD2 


dD 

dD 

az 


2  B(D) 
3  B'(D) 


B(D) 


U(Z ,  D)  =  Q(Z ,  D). 

(4.48) 


Kinetic  coefficients  LAB  are  not  needed.  Substituting 
from  (4.14)  and  (4.40),  the  two  balance  laws  (4.47) 
and  (4.48)  become,  for  the  equilibrium  downstream 
shocked  state, 


d  P  d£ 

—  =  -kP (4.49) 
dZ  dZ 

+  2 Tlf  -  -  2A£(1  -  m  -  2$) 


=  2  k% 


U  - 


r/(^)2]  -2pI(fe) 


e,x  9^D)f 

X' 


3? 


+  P(FE)Z 


3? 


(4.50) 


Terms  involving  Weyl  factor  k  result  from  dependence 
of  the  generalized  pseudo-Finsler  metric  and  corre¬ 
sponding  Cartan’s  tensor  coefficients  on  internal  state. 
For  Riemannian  geometry,  k  =  0  and  such  terms  van¬ 
ish.  Relations  (4.49)  and  (4.50)  are  two  coupled  non¬ 
linear  differential  equations  wherein  dependent  field 


variables  PA,  (FE)“,  ( FD)\ ip,  §,  and  U  all  are  ulti¬ 
mately  functions  of  independent  variable  Z.  However, 
the  downstream  state  is  assumed  to  be  spatially  homo¬ 
geneous  (in  material  coordinates)  with  regard  to  field 
variables,  meaning  that  gradients  with  respect  to  Z 
vanish  identically.  Thus,  differential  equation  (4.49) 
reduces  to  the  trivial  condition 


Q 

P  =  —  F3  =  constant  (downstream  equilibrium), 

(4.51) 


with  the  value  of  this  component  and  the  equal  and 
constant  lateral  stresses  obeying  constitutive  equations 
in  (4.39).  Since  a  homogeneous  order  parameter  field 
over  the  shocked  portion  of  201  is  in  effect,  £'  =  0 
and  =  0.  Equation  (4.50)  then  degenerates  to  the 
algebraic  equation 


-2y$-2A$(l-$)(l-2$) 

i  f  x  3 (FD)* 
=  2k'%U  -  2P\(Fe)\— - ^ 


+P(Fe)zz 


d(FD)zz 

3? 


(downstream  equilibrium). 

(4.52) 


If  compressive  deformation  6  =  F^  —  l  is  prescribed  as 
the  condition  denoting  the  intensity  of  the  shock  load¬ 
ing,  equations  (4.39),  (4.46),  and  (4.52)  can  be  solved 
simultaneously  for  the  downstream  state,  where  the 
form  of  internal  energy  in  (2.65),  (4.31),  and  (4.32) 
is  also  invoked.  In  other  words,  volume  reduction  is 
applied  incrementally,  and  then  the  balance  laws,  jump 
conditions,  and  constitutive  equations  are  solved  simul¬ 
taneously  for  the  stresses,  entropy,  internal  energy,  and 
order  parameter  in  the  material  behind  the  shock,  as 
well  as  the  particle  velocity  v  and  the  shock  speed 
D.  Since  ^(0)  =  0  from  (4.23),  the  entire  right  side 
of  (4.52)  vanishes  when  $  =  0  according  to  (4.41). 
Because  the  left  side  also  vanishes  for  $  =  0,  the  non¬ 
linear  elastic  solution  (i.e.,  no  inelastic  deformation, 
no  order  parameter  evolution)  is  clearly  always  a  solu¬ 
tion  to  (4.52).  However,  the  nonlinear  elastic  solution 
is  not  necessarily  the  only  solution,  and  it  tends  to  be 
metastable,  i.e.,  of  higher  total  energy,  than  the  alter¬ 
native  solution  that  exists  above  some  nonzero  strain 
and  for  which  $  >0.  For  the  present  application  to 
homogeneous  downstream  states,  solutions  are  sought 
in  practice  by  prescribing  a  value  of  $  =  1  as  an  initial 
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guess  at  each  applied  strain  increment  and  then  decreas¬ 
ing  §  iteratively  until  the  governing  equations  are  sat¬ 
isfied  simultaneously  to  within  a  tight  tolerance.  For 
more  general  applications  requiring  advanced  numer¬ 
ical  methods,  it  may  be  necessary  to  institute  a  per¬ 
turbation  in  §  somewhere  in  the  domain  as  an  initial 
condition  to  enable  the  inelastic  solution,  as  has  been 
invoked  in  phase  field  simulations  of  twinning  (Clayton 
and  Knap  201  la, b,  2013)  that  demonstrate  a  similar 
condition. 

The  terms  proportional  to  the  gradient  of  the  order 
parameter,  £'  =  d§/dZ,  in  (4.47)  and  (4.49)  are  ulti¬ 
mately  a  result  of  application  of  Rund’s  divergence 
theorem  to  (2.60)  for  pseudo-Finsler  space,  leading 
to  contributions  from  the  product  of  the  trace  of  Car- 
tan’s  tensor  and  the  state  variable  gradient  to  the  bal¬ 
ance  of  linear  momentum  (Clayton  2017a).  Physically, 
such  terms  can  be  interpreted  as  material  or  configu¬ 
rational  forces  induced  by  microstructure  heterogene¬ 
ity,  since  they  have  a  similar,  but  not  identical,  effect 
on  the  linear  momentum  balance  as  those  emerging 
from  Eshelby-type  forces  (e.g.,  forces  induced  by  local 
gradients  in  moduli)  in  elasticity  theory  (Marsden  and 
Hughes  1994;  Clayton  2011). 

4.3  Material  characteristics 

Boron  carbide  (B4C)  is  the  particular  material  toward 
which  the  analysis  of  Sects.  4.1  and  4.2  is  directed. 
In  its  ambient  solid  state,  boron  carbide  is  a  low  den¬ 
sity  crystalline  ceramic  of  high  hardness,  high  elastic 
stiffness,  and  low  ductility.  The  usual  crystal  structure 
is  rhombohedral.  Material  parameters  used  in  the  con¬ 
stitutive  model,  with  supporting  references,  are  listed 
in  Table  1.  If  no  reference  or  equation  is  listed  for  a 
particular  value,  the  reference  quoted  above  applies. 

Most  physical  properties  are  self-explanatory, 
though  those  in  the  final  six  rows  merit  further  dis¬ 
cussion.  The  surface  energy  is  that  corresponding  to 
fracture  on  { lOl  1}  pyramidal  planes  of  the  single  crys¬ 
tal  as  computed  via  first  principles  density  functional 
theory  (DFT)  (Beaudet  et  al.  2015).  The  regulariza¬ 
tion  length  l  and  intrinsic  surface  energy  are  specified 
to  have  magnitudes  corresponding  to  those  for  frac¬ 
ture  since  failure  accompanies  amorphization  in  exper¬ 
iments  and  since  widths  of  amorphous  zones  observed 
experimentally  are  on  the  order  of  a  nanometer  (Yan 
et  al.  2009;  Grady  2011),  of  the  same  order  as  the  frac¬ 


ture  process  zone  length.  The  value  of  /  in  Table  1  is 
computed  as  the  cohesive  fracture  process  zone  size 
over  which  the  stress  at  a  (mode  II)  crack  tip  degrades 
(Rice  1968;  Clayton  et  al.  2012;  Clayton  2017a): 

l  =  4/iT/[(l  -  v)tto2\.  (4.53) 

Here,  v  is  Poisson’s  ratio  and  a  ~  ^  is  the  theoreti- 
cal  shear  strength  of  the  crystal  (Clayton  201 1).  In  the 
intended  interpretation  of  the  theory,  as  §  — >  1  locally, 
the  material  progresses  from  glassy  sheared  state  to 
a  failed  state  at  which  inelastic  deformation  FD  satu¬ 
rates.  Subsequently,  the  material  is  expected  to  undergo 
cavitation  or  some  other  means  of  fracture.  Though  not 
implemented  in  the  present  paper,  a  reduction  in  elas¬ 
tic  coefficients  would  be  needed  to  completely  repre¬ 
sent  the  fully  failed  state,  as  in  phase  field  theories  of 
fracture  (Clayton  and  Knap  2014,  2015b).  The  surface 
energy  Y,  which  ultimately  affects  the  regularization 
length  /  through  (4.53),  is  assumed  equal  to  the  frac¬ 
ture  energy,  even  though  the  material  may  be  undergo¬ 
ing  shear  localization  rather  than  (mode  II)  fracture  for 
§  e  (0,  1).  This  assumption  is  made  in  part  because  the 
surface  energy  of  the  amorphous  bands  in  the  material 
is  not  well  known.  Thus,  the  same  energy  is  employed 
to  regularize  a  related  mechanism  (shear  localization) 
as  well  as  govern  fracture  itself  when  §  reaches  unity. 

The  Weyl  transformation  factor  k  is  determined  from 
the  ratio  of  mass  density  of  the  crystalline  phase  to  that 
of  the  glassy  phase.  This  value  is  determined  for  the 
Finsler-geometric  theory  via  consideration  of  (4.17)  at 
§  =  1 .  As  in  Clayton  (2014a),  invoked  here  is  a  4%  vol¬ 
ume  reduction  (i.e.,  mass  density  increase)  upon  struc¬ 
ture  collapse  commensurate  with  complete  amorphiza¬ 
tion  (Yan  et  al.  2009;  Taylor  2015;  An  and  Goddard 
2015a),  leading  to  exp (k /  2)  =  0.96  =>  k  =  2  In (0.96), 
consistent  with  (4.24).  The  barrier  for  phase  transfor¬ 
mation  for  a  double  well  potential  is  at  §  =0.5.  This 
barrier  is  chosen  as  the  difference  between  ground  state 
energy  of  the  most  stable  B4C  poly  type  and  the  energy 
of  segregated  elemental  phases  (boron  and  amorphous 
carbon)  associated  with  structure  collapse,  0.04  eV 
obtained  from  DFT  (Fanchini  et  al.  2006).  Finally, 
the  inelastic  shear  strain  yo  accommodated  by  amor¬ 
phous  slip  bands  at  the  onset  of  cavitation  is  obtained 
from  results  of  atomic  simulations,  specifically  molec¬ 
ular  dynamics  simulations  with  reactive  force  fields 
(RexFF)  (An  and  Goddard  2015a).  In  these  simula¬ 
tions,  simple  shearing  at  fixed  volume  on  a  pyrami- 
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Table  1  Physical 
properties  of  boron  carbide 


Property  (Units) 

Value 

Description 

References 

Cl!  (GPa) 

543 

Second-order  elastic 

constant 

Clayton  (2012a) 

C12  (GPa) 

131 

C13  (GPa) 

64 

C14  (GPa) 

-18 

C33  (GPa) 

535 

C44  (GPa) 

165 

B0  (GPa) 

237 

Bulk  modulus 

K 

4.7 

Pressure  derivative  of  bulk 
modulus 

Go  (GPa) 

236 

Shear  modulus 

c,  a  (nm) 

1.21,0.56 

Lattice  parameters 

Po  (g/cm3) 

2.51 

Mass  density 

c0  (MPa/K) 

2.41 

Specific  heat 

Clayton  and  Tonge  (2015) 

To 

1.28 

Griineisen  parameter 

Oo(K) 

295 

Ambient  temperature 

T  (J/m2) 

3.27 

Surface  energy 

Beaudet  et  al.  (2015) 

/  (nm) 

0.97 

Regularization  length 

(4.53) 

exp(&/2) 

0.96 

volume  reduction 
(amorphization) 

Yan  et  al.  (2009),  Clayton  (2014a), 
An  and  Goddard  (2015  a) 

A  (GPa) 

3.01 

Transformation  barrier 

Fanchini  et  al.  (2006) 

Yo 

1 

12 

Inelastic  shear 
accommodation 

An  and  Goddard  (2015a) 

dal  plane,  specifically  [Il01](0lTT),  resulted  in  a  total 
difference  in  shear  strain  of  ~  \  from  the  strain  at 
emergence  of  amorphous  bands  on  this  plane  at  peak 
shear  stress,  followed  by  somewhat  gradual  softening 
behavior,  and  then  finally  failure  by  cavitation  and  shear 
fracture.  Taking  the  emergence  point  in  the  Finsler  con¬ 
tinuum  representation  as  that  corresponding  to  initia¬ 
tion  of  nonzero  §,  and  taking  the  failure  point  to  corre¬ 
spond  to  §  —>  1 ,  slip  accommodation  for  a  single  plane 
would  thus  be  ^  at  saturation,  i.e.,  at  complete  transfor¬ 
mation  just  prior  to  rupture.  In  the  present  problem  of 
compression  along  [0001],  since  six  pyramidal  planes 
support  identical  amorphous  bands,  it  is  assumed  that 
degradation  resulting  from  each  plane  is  cumulative, 
leading  to  yo  =  \  •  \  =  ^ .  The  following  caveat  is 
noted,  however.  In  atomic  simulations  reported  in  An 
and  Goddard  (2015a),  steadily  increasing  shear  stresses 
reaching  local  maxima  in  excess  of  35-45  GPa,  depend¬ 
ing  on  possible  volume  relaxation,  were  observed,  with 
substantial  amorphization  not  taking  place  until  shear 
strains  in  excess  of  0.3  were  applied.  These  magnitudes 


of  shear  stress  and  shear  strain  far  exceed  those  pre¬ 
dicted  later  in  Sect.  4.4  for  shock  compression.  How¬ 
ever,  boundary  conditions  differ  substantially  for  the 
present  work  versus  (An  and  Goddard  2015a),  and  the 
shear  strengths  reported  in  the  latter  are  essentially 
upper  bounds  (i.e.,  theoretical  maxima)  since  initial 
imperfections  in  the  material  are  excluded  by  design. 
Acquisition  of  the  present  value  of  yo  motivated  from 
An  and  Goddard  (2015a)  is  perhaps  inconsistent  in  light 
of  the  above  discrepancies,  but  the  present  choice  is 
deemed  favorable  to  simply  treating  the  constant  as  an 
adjustable  parameter.  Furthermore,  the  present  choice 
is  demonstrated  later  to  enable  accurate  prediction  of 
shock  stress  versus  experimental  data. 

Importantly,  in  the  present  application  of  the  Finsler 
theory  to  shock  compression  of  boron  carbide,  param¬ 
eter  fitting  by  matching  results  of  the  model  to  exper¬ 
iments  or  third-party  simulations  is  not  necessary.  All 
values  shown  in  Table  1  are  prescribed  a  priori  and 
directly  from  experimental  data  or  atomic  calculations 
reported  in  the  literature,  though  the  aforementioned 
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Fig.  1  Boron  carbide, 
planar  shock  compression 
along  c-axis  a  axial  (shock) 
stress  and  experimental  data 
Vogler  et  al.  (2004)  b  shear 
stress  and  experimental  data 
Vogler  et  al.  (2004)  c 
Cauchy  pressure  d  internal 
energy  e  shock  velocity 
versus  particle  velocity  f 
order  parameter  f  =  D/ 1 


v/CL  1-v/V 

(e)  (f) 


caveat  regarding  yo  is  recalled.  Thus,  model  outcomes 
reported  later  in  Sect.  4.4  are  considered  fully  predic¬ 
tive. 


4.4  Solutions  and  interpretations 

Solutions  to  the  planar  shock  problem  for  boron  carbide 
are  illustrated  in  Fig.  1  which  is  analyzed  in  detail  in 
the  following  discussion.  Axial  shock  stress  P  normal¬ 
ized  by  the  bulk  modulus  Bo  is  shown  versus  compres¬ 
sion  in  Fig.  la,  where  J  =  v/V  =  F33  =  1  +  e  is  the 


ratio  of  volume  after  compression  to  the  initial  volume. 
Results  for  the  physically  realistic  case  invoking  inelas¬ 
tic  volume  reduction  in  conjunction  with  amorphiza- 
tion,  i.e.,  densification,  correspond  to  the  Weyl  scaling 
parameter  k  =  2  In (0.96)  as  listed  in  Table  1.  Results 
for  k  =  0  omit  the  density  difference  between  crystal 
and  glass  phases.  Experimental  data  correspond  to  pla¬ 
nar  impact  tests  on  the  poly  crystalline  ceramic  (Vogler 
et  al.  2004);  experimental  shock  data  for  single  crystals 
is  absent  in  the  literature.  The  isentrope  corresponding 
to  purely  elastic  uniaxial  strain  of  the  single  crystal  is 
obtained  by  setting  §  =  0  throughout  a  correspond- 
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in g  static  compression  process.  Excellent  agreement 
between  predictions  of  the  present  model  with  k  <  0 
and  experimental  data  is  obvious,  thereby  lending  con¬ 
fidence  to  the  current  theory  and  solution,  at  least  for 
this  particular  application.  The  Finsler  solutions  and 
experiment  all  suggest  a  Hugoniot  Elastic  Limit  (HEL, 
i.e.,  stress  P  at  the  initial  yield  point)  of  P  ~  18  GPa 
at  v/V  ~  0.96.  Post-yielding,  shock  stress  is  higher 
for  the  case  when  inelastic  density  changes  are  omitted 
(k  =  0)  than  when  they  are  included  (k  <  0).  Stress  P 
is  largest  for  the  isentrope,  as  expected,  since  inelastic 
deformation  by  shear  localization  and  volume  collapse 
that  would  otherwise  relieve  elastic  strain  is  excluded. 

Dynamic  shear  strength  x  normalized  by  the  shear 
modulus  Go  =  C 33  —  C13)  is  shown  versus  com¬ 
pression  in  Fig.  lb,  where  specifically 

T  =  =  -4\  (4-54) 

is  half  the  difference  between  longitudinal  and  trans¬ 
verse  components  of  Cauchy  stress  a .  Close  agreement 
between  the  present  solution  with  k  <  0  and  experi¬ 
ment  (Vogler  et  al.  2004)  is  evident.  Shear  stress  at 
the  HEL  is  x  ~  7  GPa  (Vogler  et  al.  2004);  strength 
degrades  then  increases  slightly  with  increasing  com¬ 
pression  beyond  this  point  for  the  present  model  with 
k  <  0  as  well  as  experiment.  As  is  the  case  for  axial 
stress,  shear  stress  is  highest  at  larger  compression 
for  the  isentropic  solution  (no  shear  accommodation 
by  slip/inelasticity),  with  x  for  the  Finsler  model  with 
k  =  0  falling  in  between  the  other  two  cases  at  com¬ 
pressive  strains  exceeding  that  at  the  HEL.  As  evi¬ 
denced  by  the  drop  in  shear  strength  upon  attainment 
of  the  HEL,  both  the  experimental  data  and  the  Finsler 
solutions  demonstrate  how  boron  carbide  softens  and 
fails  in  shear  under  compressive  loading.  In  the  real 
material,  at  large  compressive  deformations,  friction  at 
internal  surfaces  may  contribute  to  shear  strength;  this 
effect  is  omitted  in  the  present  application  of  the  theory. 

A  limitation  of  the  present  shock  solutions  is  the 
predicted  deviation  from  realistic  behavior  at  volumet¬ 
ric  compressions  exceeding  20%,  for  which  an  upturn 
in  shear  stress  is  evident  in  Fig.  lb.  This  upturn  is  a 
result  of  the  contribution  from  internal  state  §  to  (4.22) 
approaching  saturation,  such  that  inelastic  shear  defor¬ 
mation  is  no  longer  able  to  sufficiently  offset  total  devi- 
atoric  deformation,  leading  to  an  increase  in  elastic 
strain  and  shear  stress.  The  curvature  of  the  curve  of 


total  axial  stress  P  in  Fig.  1  a  likewise  becomes  too  great 
at  large  compression.  The  problem  could  be  remedied 
by  permitting  elastic  stiffness  coefficients  (e.g.,  shear 
moduli  but  not  compressive  bulk  modulus)  to  degrade 
with  increasing  §,  as  in  phase  field  fracture  models 
(Clayton  and  Knap  2014,  2015a),  above  some  thresh¬ 
old.  Such  model  additions  should  be  enabled  in  future 
work  for  more  general  materials  failure  problems. 

Figure  lc  reports  the  Cauchy  pressure  p  normalized 
by  the  ambient  bulk  modulus,  computed  via 

P  =  +  al  +  ct3)  =  +  CT3)  =  P  ~2t, 

(4.55) 

recalling  that  P  =  —  a33.  The  same  three  cases  are 
addressed:  the  Finsler  model  with  no  phase  densifica- 
tion  (k  =  0),  the  Finsler  model  with  realistic  inelastic 
densification  [k  =  2  In (0.96)],  and  isentropic  compres¬ 
sion  (§  =0).  Complementary  predictions  for  the  inter¬ 
nal  energy  per  unit  reference  volume  U  are  shown  in 
Fig.  Id.  Results  of  the  three  cases  coincide  for  compres¬ 
sive  loading  below  the  HEL,  i.e.,  for  v/V  >  0.96.  For 
larger  compression,  pressure  is  largest  for  the  Finsler 
result  with  k  =  0,  followed  by  the  isentrope,  and  then 
smallest  for  the  Finser  result  with  k  <  0.  Comparing 
Fig.  la  with  c,  the  shock  stress  P  only  slightly  exceeds 
the  Cauchy  pressure  p  for  compressive  volume  reduc¬ 
tions  in  the  range  0.95  >  v/V  >  0.80  where  shear 
stress  is  relatively  small  (Fig.  lb).  Ordering  of  internal 
energy  is  interchanged  between  the  cases  of  isentrope 
and  Finsler  result  with  k  —  0.  In  particular,  results  for 
internal  energy  demonstrate  that  the  Finsler  solutions 
are  stable,  i.e.,  of  lower  total  internal  energy,  relative 
to  the  isentropic  elastic  solution  that  has  the  largest 
energy.  This  is  an  important  finding  since  the  isentropic 
solution  is  always  a  (possibly  non-unique)  solution  to 
the  static  governing  equations  in  the  present  constitu¬ 
tive  model.  Pressure  is  larger  in  the  Finsler  result  for 
k  —  0  than  that  of  the  isentropic  solution  because  of 
thermoelastic  coupling  and  substantial  entropy  produc¬ 
tion.  The  increase  in  pressure  from  entropy  production 
is  more  than  offset  by  a  decrease  in  pressure  with  den¬ 
sification  for  the  Finsler  result  with  k  =  2  In (0.96). 
Normalized  temperatures  6/60  at  20%  compression  for 
the  shock  solutions  are  4.0  and  4.1  for  k  =  0  and 
k  =  2  In (0.96),  compared  to  1 .3  for  the  isentropic  solu¬ 
tion.  The  temperature  rise  and  pressure  for  the  shocked 
material  in  the  absence  of  inelasticity  (i.e.,  for  a  non¬ 
linear  elastic  shock,  physically  valid  only  for  shock 
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stresses  up  to  the  HEL)  only  slightly  exceed  those  of 
the  isentrope;  this  has  been  demonstrated  elsewhere  in 
calculations  for  a  number  of  brittle  ceramics  and  min¬ 
erals  invoking  nonlinear  logarithmic  thermoelasticity 
(Clayton  2014c,  b;  Clayton  and  Tonge  2015).  A  con¬ 
tribution  to  tensile  pressure  in  amorphous  bands  lead¬ 
ing  to  cavitation  has  been  noted  in  prior  atomic  (An 
and  Goddard  2015a)  and  continuum  (Clayton  2016b) 
model  results  for  shear  loading. 

Shock  velocity  D  is  shown  versus  particle  velocity 
v  in  the  compressed  state  in  Fig.  le,  with  both  veloci¬ 
ties  normalized  by  the  longitudinal  elastic  wave  speed 
Cl  =  yC33/po-  The  same  three  cases  are  considered, 
with  the  isentropic  solution  giving  a  steadily  increas¬ 
ing  shock  speed  with  increasing  particle  velocity  corre¬ 
sponding  to  increasing  shock  stress  or  decreasing  vol¬ 
ume.  This  increasing  isentropic  shock  speed  is  a  result 
of  elastic  nonlinearity,  specifically  the  increasing  tan¬ 
gent  longitudinal  modulus  with  increasing  elastic  com¬ 
pression.  When  the  shock  stress  exceeds  the  HEL  for 
the  other  two  cases,  a  drop  in  shock  velocity  occurs 
commensurate  with  a  loss  of  tangent  stiffness  and  shear 
strength  (Fig.  lb).  The  largest  drop  in  shock  veloc¬ 
ity  occurs  for  the  realistic  model  with  k  <  0.  Shock 
velocity  reaches  a  local  minimum  then  increases  with 
increasing  compressive  strain  or  particle  velocity.  The 
same  trends  in  shock  velocity  versus  particle  velocity 
are  evident  in  experimental  data  of  Vogler  et  al.  (2004) 
and  in  many,  but  not  all,  experimental  and  quantum 
molecular  results  reported  in  Taylor  (2015). 

Finally,  the  order  parameter  physically  connected  to 
shear  accommodation  by  amorphous  bands  is  shown 
versus  compressive  strain  in  Fig.  If  for  the  Finsler 
model  with  k  =  0  and  k  =  2  In (0.96).  At  small 
compressive  strains  v/V  >  0.97,  the  order  parame¬ 
ter  §  =0  and  the  response  is  isentropic  and  elastic.  A 
rapid  increase  occurs  for  each  case  around  the  HEL, 
followed  by  a  more  gradual  increase  at  larger  com¬ 
pressive  strains  or  higher  shock  stresses.  The  magni¬ 
tude  of  §  for  the  case  with  densification  upon  phase 
transformation  slightly  exceeds  that  for  the  case  with 
k  =  0  since  the  magnitude  of  the  elastic  driving  force 
3  W/d%  promoting  evolution  of  §  is  greater  under  com¬ 
pressive  loading  when  volume  reduction  due  to  amor- 
phization  is  included  in  the  constitutive  model.  In  other 
words,  volume  reduction  due  to  phase  transformation 
accommodates  strain  that  otherwise  would  be  accom¬ 
modated  elastically  leading  to  an  increase  in  pressure 
and  volumetric  strain  energy  density.  Notice  that  in 


each  case,  the  order  parameter  increases  towards,  but 
never  attains,  a  value  of  unity  with  increasing  applied 
strain.  The  reason  for  this  arises  from  a  nonzero  ther¬ 
modynamic  force  (2)  contribution  from  the  quadratic 
term  in  §  entering  the  energy  functional,  leading  to  the 
first  term  on  the  left  side  of  (4.52)  that  does  not  vanish 
at  §  =  1 .  Such  an  effect  is  also  present  in  phase  field 
models  of  fracture  that  employ  the  standard  quadratic 
contribution  from  an  order  parameter  to  the  energy  den¬ 
sity  (Clayton  and  Knap  2015a;  Borden  et  al.  2012). 
Homogeneous  solutions  for  quasi-static  fracture  prob¬ 
lems  demonstrate  order  parameter  values  that  asymp¬ 
totically  approach  unity  with  increasing  strain,  with  the 
rate  of  approach  dependent  upon  normalized  regular¬ 
ization  length  (Clayton  and  Knap  2015a). 

Recent  work  involving  either  quantum  mechanics 
(DFT)  or  molecular  dynamics  simulations  (ReaxFF) 
(An  and  Goddard  2015a,  b;  Tang  et  al.  2015)  has  pro¬ 
vided  insight  into  effects  of  composition  and  load¬ 
ing  direction  (with  respect  to  the  crystal  structure) 
on  resistance  to  shear  localization,  subsequent  cavita¬ 
tion,  and  failure  of  boron-based  ceramic  crystals  under 
(simple)  shear  loading  with  periodic  boundary  condi¬ 
tions.  In  particular,  computationally  designed  atomic 
structures  featuring  layers  of  B4C  and  boron  subox¬ 
ide  (B60)  have  demonstrated  increased  overall  ductil¬ 
ity  relative  to  more  brittle  pure  boron  carbide  crystals 
(Tang  et  al.  2015).  It  is  suggested  here  that  improve¬ 
ments  in  fundamental  physical  properties  such  as  sur¬ 
face  energy  T,  glassy  phase  characteristics  (e.g.,  mass 
density  change  quantified  by  k  and  energy  change  quan¬ 
tified  by  A),  and  shear  slip  accommodation  by  amor¬ 
phous  banding  quantified  by  yo  may  be  possible  via  tai¬ 
loring  of  boron-based  ceramic  compositions  and  struc¬ 
tures.  Conversely,  as  indicated  by  quantum  mechani¬ 
cal  results  (Taylor  2015),  physical  properties  such  as 
resistance  to  structure  collapse  may  be  deleteriously 
affected  via  changes  in  atomic  structure,  for  example 
lattice  site  vacancies. 

Predictions  of  the  present  continuum  model  of  shock 
compression  demonstrating  effects  of  variations  of  the 
aforementioned  fundamental  properties  T,  A,  and  yo 
are  shown  in  Fig.  2.  Specifically,  Fig.  2a  shows  shock 
stress  P ,  the  fundamental  mechanical  force  resisting 
uniaxial  shock  compression,  while  Fig.  2b  shows  shear 
stress  (i.e.,  strength)  r ,  which  is  thought  to  be  a  primary 
indicator  of  resistance  of  a  ceramic  material  to  ballis¬ 
tic  penetration,  i.e.,  penetration  by  an  impacting  pro¬ 
jectile  (Bourne  2008;  Clayton  2015b).  Results  labeled 
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Fig.  2  Effects  of  surface 
energy  T,  transformation 
barrier  A,  and  shear  band 
strain  yo  on  a  axial  (shock) 
stress  b  shear  stress 


1-v/V  l-v/V 


(a) 


(b) 


as  “nominal”  correspond  to  the  most  physically  real¬ 
istic  model  considered  here  with  k  =  2  In (0.96),  as 
discussed  already  in  the  context  of  Fig.  1.  Recall  that 
effects  of  k  have  already  been  reported  in  Fig.  1 ,  demon¬ 
strating  for  k  =  0  an  increase  in  shock  stress  and  shear 
strength  relative  to  the  nominal  case.  General  increases 
in  shock  stress  and,  more  evidently,  shear  strength  are 
apparent  with  increases  in  fracture  energy  and  trans¬ 
formation  energy,  and  with  decreases  in  yo.  The  total 
driving  force  for  transformation  decreases  commensu- 
rately  with  such  property  variations,  thereby  inhibiting 
the  trajectory  of  the  shocked  substance  towards  a  fully 
amorphous  and  subsequently  ruptured  state.  Notewor¬ 
thy  is  the  dramatic  increase  in  strength  with  decreasing 
yo,  clarifying  remarks  from  Sect.  4.3  that  this  parame¬ 
ter  is  not  a  measure  of  ductility  of  the  crystal,  but  rather 
is  better  labeled  as  a  measure  of  the  shear  damage  asso¬ 
ciated  with  amorphous  band  formation  and  multiplica¬ 
tion,  somewhat  analogous  to  the  overall  inelastic  shear 
strain  accommodated  by  crack  opening  displacement 
for  mode  II  fracture.  The  greatest  potential  increases 
to  strength  correspond  to  halving  of  this  parameter 
and  doubling  of  surface  energy  T.  Less  improvement 
results  from  doubling  of  transformation  (double  well 
potential)  barrier  A,  which  affords  almost  no  strength 
increase  at  very  large  compressions  relative  to  the  nom¬ 
inal  case. 

The  example  considered  herein  has  imposed  homo¬ 
geneous  deformation  and  internal  state  behind  the 
shock  front,  with  total  deformation  one-dimensional. 
This  example  is  deemed  a  requisite  early  step  towards 
understanding  and  eventually  validating  a  novel  and 
mathematically  sophisticated  theory.  Future  work  will 
consider  multi-dimensional  problems  involving  hetero¬ 
geneous  deformation  and  spatial  gradients  of  state  vari¬ 
ables,  including  explicit  localization  phenomena.  More 


advanced  numerical  methods  will  be  needed  to  solve 
such  problems  whose  solutions  should  provide  addi¬ 
tional  insight  into  material  performance,  i.e.,  resistance 
to  failure. 


5  Conclusions 

A  theory  of  Finsler-geometric  continuum  mechanics 
developed  recently  by  the  author  using  variational 
principles  has  been  extended  to  address  dynamics — 
material  inertia,  order  parameter  evolution,  internal 
energy  conservation — as  well  as  temperature  change 
and  entropy  production  for  the  adiabatic  case.  Jump 
conditions  pertinent  to  exchange  of  mass,  linear 
momentum,  and  energy  across  a  shock  front  moving 
at  steady  speed  have  been  derived.  The  theory  has  been 
invoked  to  describe  stress-induced  phase  transitions 
and  shear  inelasticity  in  boron  carbide  subjected  to  pla¬ 
nar  impact-type  loading.  All  parameters  in  the  model 
are  obtained  from  fundamental  experiments  or  results 
of  atomic  simulations  from  the  literature,  without  fur¬ 
ther  calibration.  Predictions  of  the  model  for  shock 
stress,  shear  strength,  and  shock  characteristics  agree 
closely  with  experimental  data  when  the  conformal 
transformation  accounting  for  densification  upon  phase 
change  is  included.  In  order  to  guide  efforts  towards 
design  of  new  compositions  of  boron-based  ceramics, 
parameter  studies  predicting  effects  of  variations  in  sur¬ 
face  energy,  transformation  energy,  and  shear  accom¬ 
modation  have  been  reported.  Results  suggest  that,  in 
increasing  order  of  improved  dynamic  shear  strength, 
the  following  structural  property  changes  should  be 
sought:  increases  in  the  energy  barrier  for  amorphiza- 
tion,  decreases  in  the  density  of  the  glassy  phase, 


Springer 


Finsler-geometric  continuum  dynamics  and  shock  compression 


77 


increases  in  the  surface  energy,  and  decreases  in  the 
post-peak  shear  localization  strain. 
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